Does philopatry pay? Meta-analysis reveals little empirical support for dispersal being costly.

Supplmentary Material

Things to do for April

1. Rerun all the models
2. Put figures at the end
3. Add some text so the reader can understand
4. Go through code and add more annotations
5. Ask questions about things you do not understand

1 Setting up

1.1 Load packages

pacman::p_load(tidyverse, # tidy family and related pacakges below
               kableExtra, # nice tables
               MuMIn,  # multi-model inference
               emmeans, # post-hoc tests
               gridExtra, # may not use this
               pander,   # nice tables
               metafor,  # package for meta-analysis
               ape,      # phylogenetic analysis
               MuMIn,  # multi-model inference
               patchwork,   # putting ggplots together - you need to install via devtool
               here,         # making reading files easy
               orchaRd, # plotting orchard plots
               matrixcalc # matrix calculations
               # more too add
)

eval(metafor:::.MuMIn)

1.2 Load data: paper and tree data

#the dataset with effect sizes
dat <- read.csv(here("data", "clean_data.csv"))

#dim(dat)
#head(dat)

#the phylogenetic tree
tree <- read.tree(here("data", "species_tree.tre"))

1.3 Calculate effect sizes and sampling variances using custom functions

# function to calculate effect sizes
# Zr - correlation
# there is always n

effect_size <- function(m1, m2, sd1, sd2, n1, n2, n, # 12 arguments
                        est , se, p_val, direction, method){

  if(method == "mean_method"){
  
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s_pool <- sqrt( ((n1-1)*sd1^2 + (n2-1)*sd2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r    
    
  }else if(method == "count_method"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s1 <- sqrt(m1)
    s2 <- sqrt(m2)
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r     
    
  }else if(method == "percent_method1"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    m1 <- m1/100
    m2 <- m2/100
    s1 <- 1/sqrt(8)
    s2 <- 1/sqrt(8)
    m1 <- asin(sqrt(m1/100))
    m2 <- asin(sqrt(m2/100))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "percent_method2"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    m1 <- m1/100
    m2 <- m2/100
    sd1 <- sd1/100
    sd2 <- sd2/100
    s1 <- 1/sqrt(sd1^2/(4*m1*(1-m1)))
    s2 <- 1/sqrt(sd2^2/(4*m2*(1-m2)))
    m1 <- asin(sqrt(m1/100))
    m2 <- asin(sqrt(m2/100))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "proportion_method1"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s1 <- 1/sqrt(8)
    s2 <- 1/sqrt(8)
    m1 <- asin(sqrt(m1))
    m2 <- asin(sqrt(m2))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "proportion_method2"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s1 <- 1/sqrt(sd1^2/(4*m1*(1-m1)))
    s2 <- 1/sqrt(sd2^2/(4*m2*(1-m2)))
    m1 <- asin(sqrt(m1/100))
    m2 <- asin(sqrt(m2/100))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "t_method1"){
    

    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- est/sqrt(est^2 + n - 2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
    
  }else if(method == "t_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- est/sqrt(est^2 + n - 2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
    
  }else if(method == "F_method1"){
    
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- sqrt(est)/sqrt(est + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b = r_b*(direction)
    r <- r_b
  
  }else if(method == "F_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- sqrt(est)/sqrt(est + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b = r_b*(direction)
    r <- r_b
    
  }else if(method == "p_method1"){
    
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- qt(1 - p_val, n - 2)
    r_pb <- t/sqrt(t^2 + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b <- r_b*(direction)
    r <- r_b
    
  }else if(method == "p_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- qt(1 - p_val, n - 2)
    r_pb <- t/sqrt(t^2 + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b <- r_b*(direction)
    r <- r_b
    
  }else if(method == "correlation_method1"){
    
    r <- est
    
  }else if(method == "correlation_method2"){
    
    r <- 2*sin((pi/6)*est)
    
  }else if(method == "estimate_method1"){
    
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- est/se
    r_pb <- t/sqrt(t^2+ n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
  
  }else if(method == "estimate_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- est/se
    r_pb <- t/sqrt(t^2+ n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
  
  } 
  
    if(r >= 1){
    # if over 1, we use 0.99
    Zr <- atanh(0.99)

    }else if(r <= -1){

    Zr <- atanh(-0.99) # r = correlation

    } else {

    Zr <- atanh(r) # r = correlation

    }
  
  VZr <- 1 /(n - 3)
  
  data.frame(ri = r, yi = Zr , vi = VZr)
  
}


##########

#' Title: Contrast name generator
#'
#' @param name: a vector of character strings
cont_gen <- function(name) {
  combination <- combn(name, 2)
  name_dat <- t(combination)
  names <- paste(name_dat[, 1], name_dat[, 2], sep = "-")
  return(names)
}

#' @title get_pred1: intercept-less model
#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)
#' @param model: rma.mv object 
#' @param mod: the name of a moderator 
get_pred1 <- function (model, mod = " ") {
  name <- firstup(as.character(stringr::str_replace(row.names(model$beta), mod, "")))
  len <- length(name)
  
   if (len != 1) {
        newdata <- diag(len)
        pred <- metafor::predict.rma(model, 
                                     newmods = newdata,
                                     tau2.levels = 1:len)
    }
    else {
        pred <- metafor::predict.rma(model)
  }
  estimate <- pred$pred
  lowerCL <- pred$ci.lb
  upperCL <- pred$ci.ub 
  lowerPR <- pred$cr.lb
  upperPR <- pred$cr.ub 
  
  table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = model$pval,
                  lowerPR = lowerPR, upperPR = upperPR)
}

#' @title get_pred2: normal model
#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)
#' @param model: rma.mv object 
#' @param mod: the name of a moderator 
get_pred2 <- function (model, mod = " ") {
  name <- as.factor(str_replace(row.names(model$beta), 
                                paste0("relevel", "\\(", mod,", ref = name","\\)"),""))
  len <- length(name)
  
  if(len != 1){
  newdata <- diag(len)
  pred <- predict.rma(model, intercept = FALSE, newmods = newdata[ ,-1])
  }
  else {
    pred <- predict.rma(model)
  }
  estimate <- pred$pred
  lowerCL <- pred$ci.lb
  upperCL <- pred$ci.ub 
  lowerPR <- pred$cr.lb
  upperPR <- pred$cr.ub 
  
  table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = model$pval,
                  lowerPR = lowerPR, upperPR = upperPR)
}

#' @title mr_results
#' @description Function to put results of meta-regression and its contrasts
#' @param res1: data frame 1
#' @param res1: data frame 2
mr_results <- function(res1, res2) {
  restuls <-tibble(
    `Fixed effect` = c(as.character(res1$name), cont_gen(res1$name)),
    Estimate = c(res1$estimate, res2$estimate),
    `Lower CI [0.025]` = c(res1$lowerCL, res2$lowerCL),
    `Upper CI  [0.975]` = c(res1$upperCL, res2$upperCL),
    `P value` = c(res1$pval, res2$pval),
    `Lower PI [0.025]` = c(res1$lowerPR, res2$lowerPR),
    `Upper PI  [0.975]` = c(res1$upperPR, res2$upperPR),
  )
}


#' @title all_models
#' @description Function to take all possible models and get their results
#' @param model: intercept-less model
#' @param mod: the name of a moderator 

all_models <- function(model, mod = " ", type = "normal") {
  
  # getting the level names out
  level_names <- levels(factor(model$data[[mod]]))
  dat2 <- model$data
  mod <- mod

  # creating a function to run models
  run_rma1 <- function(name) {
      
    # variance covarinace matrix for sampling errors
    VCV1 <- vcalc(vi = dat$vi, cluster = dat$shared_group, rho = 0.5)
    VCV1 <- nearPD(VCV1)$mat
    
      
    rma.mv(yi, V = VCV1,
                   mods = ~relevel(dat2[[mod]], ref = name),
                     random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned),
                     data = dat2,
                     test = "t",
                     sparse = TRUE,
                     control = list(optimizer = "Nelder-Mead"))
   }

    run_rma2 <- function(name) {
    
            # variance covarinace matrix for sampling errors
            VCVa <- vcalc(vi = dat2$vi_abs, cluster = dat$shared_group, rho = 0.5)
            VCVa <- nearPD(VCVa)$mat
               
               rma.mv(abs_yi2, V = VCVa,
               mods = ~relevel(dat2[[mod]], ref = name),
                      random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned),
                     data = dat2,
                     test = "t",
                     sparse = TRUE,
                     control = list(optimizer = "Nelder-Mead"))
    }
    
    run_rma3 <- function(name) {
      
              # variance covarinace matrix for sampling errors
                VCV1 <- vcalc(vi = dat$vi, cluster = dat$shared_group, rho = 0.5)
                VCV1 <- nearPD(VCV1)$mat
               
                rma.mv(yi, V = VCV1,
                   mods = ~relevel(dat2[[mod]], ref = name),
                     random = list( # hetero in relation to mod
                            formula(paste("~", mod, "| effectID")),
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned),
                     struct = "DIAG",
                     data = dat2,
                     test = "t",
                     sparse = TRUE,
                     control = list(optimizer = "Nelder-Mead"))
   }


# results of meta-regression including all contrast results; taking the last level out ([-length(level_names)])
# use different specifications for aboslute and hetero models 
if (type == "normal"){

    model_all <- purrr::map(level_names[-length(level_names)], run_rma1)

  }else if(type == "absolute"){
    model_all <- purrr::map(level_names[-length(level_names)], run_rma2)
    
  }else if(type == "hetero"){
    model_all <- purrr::map(level_names[-length(level_names)], run_rma3)
    
  }
  
  # getting estimates from intercept-less models (means for all the groups)
  res1 <- get_pred1(model, mod = mod)
  
  # getting estiamtes from all contrast models
  res2_pre <- purrr::map(model_all, ~ get_pred2(.x, mod = mod))
  
  # a list of the numbers to take out unnecessary contrasts
  contra_list <- Map(seq, from=1, to=1:(length(level_names) - 1))
  res2 <- purrr::map2_dfr(res2_pre, contra_list, ~.x[-(.y), ]) 
  # creating a table
  res_tab <- mr_results(res1, res2) %>% 
  kable("html",  digits = 3) %>%
  kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%")
  
  # results
  res_tab

}

# absolute value analyses

# folded mean
folded_mu <-function(mean, variance){
  mu <- mean
  sigma <- sqrt(variance)
  fold_mu <- sigma*sqrt(2/pi)*exp((-mu^2)/(2*sigma^2)) + mu*(1 - 2*pnorm(-mu/sigma))
  fold_mu
} 

# folded variance
folded_v <-function(mean, variance){
  mu <- mean
  sigma <- sqrt(variance)
  fold_mu <- sigma*sqrt(2/pi)*exp((-mu^2)/(2*sigma^2)) + mu*(1 - 2*pnorm(-mu/sigma))
  fold_se <- sqrt(mu^2 + sigma^2 - fold_mu^2)
  # adding se to make bigger mean
  fold_v <-fold_se^2
  fold_v
} 

1.4 Preparing final dataset

2 Intercept meta-analytic model

2.1 Phylogenetic multilevel models

With all five random effects

# takes a while to run
mod <- rma.mv(yi = yi, 
              V = VCV,
              mod = ~ 1,
              data = dat,
              random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned,
                            ~ 1 | phylogeny),
              R= list(phylogeny = cor_tree),
              test = "t",
              sparse = TRUE)

# saving the runs
saveRDS(mod, here("Rdata", "mod.RDS"))

# getting the saved model
mod <- readRDS(here("Rdata", "mod.RDS"))

summary(mod)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-279.0913   558.1826   568.1826   590.7488   568.2724   

Variance Components:

            estim    sqrt  nlvls  fixed           factor    R 
sigma^2.1  0.1105  0.3324    675     no         effectID   no 
sigma^2.2  0.0110  0.1050    202     no          paperID   no 
sigma^2.3  0.0287  0.1694    146     no  species_cleaned   no 
sigma^2.4  0.0000  0.0000    146     no        phylogeny  yes 

Test for Heterogeneity:
Q(df = 674) = 390155258.0414, p-val < .0001

Model Results:

estimate      se     tval   df    pval    ci.lb   ci.ub    
 -0.0301  0.0303  -0.9947  674  0.3202  -0.0896  0.0294    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Phylogeney accounts for nothing (0.00%), so we can take it out. 
round(i2_ml(mod),2)   

          I2_Total        I2_effectID         I2_paperID I2_species_cleaned 
             97.74              71.90               7.17              18.67 
      I2_phylogeny 
              0.00 

pred_mod <- predict.rma(mod)

  estimate <- pred_mod$pred
  lowerCL <- pred_mod$ci.lb
  upperCL <- pred_mod$ci.ub 
  lowerPR <- pred_mod$cr.lb
  upperPR <- pred_mod$cr.ub 
  
  table_mod <- tibble("Fixed effect" = "intercept", estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = mod$pval,
                  lowerPR = lowerPR, upperPR = upperPR)

# creating a table
tab_mod <- table_mod %>%
  kable("html",  digits = 3) %>%
  kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%")

# saving tab_mod as RDS
saveRDS(tab_mod, here("Rdata", "tab_mod.RDS"))

Table S1. Mean estimate with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval) for intercept meta-analytic model with all five random effects.

tab_mod <-readRDS(here("Rdata", "tab_mod.RDS"))

tab_mod

Fixed effect	estimate	lowerCL	upperCL	pval	lowerPR	upperPR
intercept	-0.03	-0.09	0.029	0.32	-0.793	0.733

Reduced multilevel model

We remove phylogeny as a random effect from our meta-analytic model because phylogeny accounts for nothing (0%)

mod1 <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)

summary(mod1)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.9876   557.9752   565.9752   584.0636   566.0345   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1095  0.3309    681     no         effectID 
sigma^2.2  0.0110  0.1049    203     no          paperID 
sigma^2.3  0.0288  0.1698    147     no  species_cleaned 

Test for Heterogeneity:
Q(df = 680) = 390817153.0648, p-val < .0001

Model Results:

estimate      se     tval   df    pval    ci.lb   ci.ub    
 -0.0293  0.0302  -0.9690  680  0.3329  -0.0885  0.0300    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

round(i2_ml(mod1),2)   

          I2_Total        I2_effectID         I2_paperID I2_species_cleaned 
             97.71              71.65               7.19              18.86 

reduced_intercept<-orchard_plot(mod1, xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

pred_mod1 <- predict.rma(mod1)

  estimate <- pred_mod1$pred
  lowerCL <- pred_mod1$ci.lb
  upperCL <- pred_mod1$ci.ub 
  lowerPR <- pred_mod1$cr.lb
  upperPR <- pred_mod1$cr.ub 
  
  table_mod1 <- tibble("Fixed effect" = "intercept", estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = mod1$pval,
                  lowerPR = lowerPR, upperPR = upperPR)

# creating a table
tab_mod1 <- table_mod1 %>%
  kable("html",  digits = 3) %>%
  kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%")

# saving tab_mod1 as RDS
saveRDS(tab_mod1, here("Rdata", "tab_mod1.RDS"))

Table S2. Mean estimate with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval) for reduced intercept meta-analytic model without phylogeny as a random effect.

tab_mod1 <-readRDS(here("Rdata", "tab_mod1.RDS"))

tab_mod1

Fixed effect	estimate	lowerCL	upperCL	pval	lowerPR	upperPR
intercept	-0.029	-0.089	0.03	0.333	-0.79	0.732

3 Uni-moderator meta-regression models

For each of our moderators, we run a uni-moderator meta-regression model

Species class

mod_class <- rma.mv(yi = yi, V = VCV,
               mod = ~ species_class - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_class)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-276.6220   553.2440   571.2440   611.8765   571.5147   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1081  0.3287    681     no         effectID 
sigma^2.2  0.0080  0.0895    203     no          paperID 
sigma^2.3  0.0397  0.1992    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 675) = 389977293.4618, p-val < .0001

Test of Moderators (coefficients 1:6):
F(df1 = 6, df2 = 675) = 0.5443, p-val = 0.7746

Model Results:

                             estimate      se     tval   df    pval    ci.lb 
species_classActinopterygii    0.0016  0.1753   0.0093  675  0.9926  -0.3425 
species_classArachnida         0.0398  0.2898   0.1373  675  0.8909  -0.5293 
species_classAves             -0.0276  0.0375  -0.7360  675  0.4620  -0.1011 
species_classInsecta           0.0948  0.1506   0.6295  675  0.5293  -0.2009 
species_classMammalia         -0.0667  0.0463  -1.4409  675  0.1501  -0.1576 
species_classReptilia          0.0681  0.1453   0.4689  675  0.6393  -0.2172 
                              ci.ub    
species_classActinopterygii  0.3458    
species_classArachnida       0.6088    
species_classAves            0.0460    
species_classInsecta         0.3904    
species_classMammalia        0.0242    
species_classReptilia        0.3535    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_class)

   R2_marginal R2_conditional 
   0.006139514    0.310349425 

orchard_plot(mod_class, mod = "species_class", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S2. Mean fitness effect of different taxonomic classes. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (blue circles) scaled by precision (1/SE).

tab_class <- all_models(mod_class,  mod = "species_class", type = "normal")

# saving tab_class as RDS
saveRDS(tab_class, here("Rdata", "tab_class.RDS"))

Table S3. Mean estimates of different taxonomic classes with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

tab_class <-readRDS(here("Rdata", "tab_class.RDS"))

tab_class

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
Actinopterygii	0.001	-0.344	0.346	0.996	-0.850	0.851
Arachnida	0.039	-0.531	0.610	0.892	-0.925	1.004
Aves	-0.027	-0.101	0.046	0.464	-0.808	0.753
Insecta	0.094	-0.202	0.391	0.533	-0.738	0.926
Mammalia	-0.071	-0.163	0.022	0.133	-0.853	0.712
Reptilia	0.067	-0.219	0.353	0.645	-0.761	0.895
Actinopterygii-Arachnida	0.038	-0.626	0.703	0.910	-0.984	1.061
Actinopterygii-Aves	-0.028	-0.378	0.321	0.873	-0.881	0.824
Actinopterygii-Insecta	0.093	-0.358	0.545	0.685	-0.806	0.993
Actinopterygii-Mammalia	-0.071	-0.423	0.280	0.690	-0.925	0.782
Actinopterygii-Reptilia	0.066	-0.376	0.508	0.769	-0.828	0.961
Arachnida-Aves	-0.067	-0.640	0.507	0.819	-1.033	0.899
Arachnida-Insecta	0.055	-0.586	0.696	0.866	-0.953	1.063
Arachnida-Mammalia	-0.110	-0.685	0.465	0.708	-1.077	0.857
Arachnida-Reptilia	0.028	-0.607	0.663	0.931	-0.976	1.032
Aves-Insecta	0.122	-0.180	0.424	0.429	-0.712	0.956
Aves-Mammalia	-0.043	-0.147	0.061	0.417	-0.827	0.741
Aves-Reptilia	0.095	-0.194	0.384	0.520	-0.735	0.924
Insecta-Mammalia	-0.165	-0.470	0.140	0.289	-1.000	0.670
Insecta-Reptilia	-0.027	-0.434	0.380	0.896	-0.905	0.850
Mammalia-Reptilia	0.138	-0.152	0.428	0.352	-0.692	0.967

Study type

mod_type <- rma.mv(yi = yi,  
               V = VCV,
               mod = ~ study_type - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_type)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.5770   557.1540   567.1540   589.7571   567.2431   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1093  0.3306    681     no         effectID 
sigma^2.2  0.0101  0.1007    203     no          paperID 
sigma^2.3  0.0312  0.1766    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 679) = 390810607.4481, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 679) = 0.6289, p-val = 0.5335

Model Results:

                        estimate      se     tval   df    pval    ci.lb   ci.ub 
study_typenatural        -0.0273  0.0310  -0.8815  679  0.3784  -0.0882  0.0336 
study_typesemi-natural   -0.0618  0.0683  -0.9058  679  0.3654  -0.1958  0.0722 
                          
study_typenatural         
study_typesemi-natural    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_type)

   R2_marginal R2_conditional 
  0.0007287837   0.2749149432 

orchard_plot(mod_type, mod = "study_type", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S3. Mean fitness effect of natural and semi-natural study types. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (blue circles) scaled by precision (1/SE).

tab_type <- all_models(mod_type,  mod = "study_type", type = "normal")

# saving tab_type as RDS
saveRDS(tab_type, here("Rdata", "tab_type.RDS"))

Table S4. Mean estimates of natural and semi-natural study types with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

tab_type <-readRDS(here("Rdata", "tab_type.RDS"))

tab_type

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
Natural	-0.028	-0.089	0.033	0.364	-0.795	0.738
Semi-natural	-0.062	-0.196	0.072	0.366	-0.838	0.714
Natural-Semi-natural	-0.034	-0.165	0.098	0.616	-0.809	0.742

Study design

mod_design <- rma.mv(yi = yi, V = VCV,
               mod = ~ study_design - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_design)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.7316   557.4631   567.4631   590.0662   567.5523   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1094  0.3308    681     no         effectID 
sigma^2.2  0.0118  0.1087    203     no          paperID 
sigma^2.3  0.0285  0.1687    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 679) = 390773683.0248, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 679) = 0.5046, p-val = 0.6040

Model Results:

                      estimate      se     tval   df    pval    ci.lb   ci.ub 
study_designcontrast   -0.0264  0.0318  -0.8307  679  0.4064  -0.0889  0.0361 
study_designgradient   -0.0389  0.0461  -0.8434  679  0.3993  -0.1294  0.0516 
                        
study_designcontrast    
study_designgradient    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_design)

   R2_marginal R2_conditional 
  0.0001882001   0.2693179740 

orchard_plot(mod_design, mod = "study_design", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S4. Mean fitness effect of gradient and contrast study designs. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (blue circles) scaled by precision (1/SE).

tab_design <- all_models(mod_design,  mod = "study_design", type = "normal")

# saving tab_design as RDS
saveRDS(tab_design, here("Rdata", "tab_design.RDS"))

Table S5. Mean estimates of gradient and contrast study designs with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

tab_design <-readRDS(here("Rdata", "tab_design.RDS"))

tab_design

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
Contrast	0.007	-0.033	0.047	0.735	-0.575	0.589
Gradient	-0.039	-0.106	0.028	0.257	-0.624	0.546
Contrast-Gradient	-0.025	-0.107	0.057	0.553	-0.758	0.709

Dispersal type

mod_disp <- rma.mv(yi = yi, 
               V = VCV,
               mod = ~ dispersal_type - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_disp)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.1612   556.3224   568.3224   595.4373   568.4476   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1097  0.3312    681     no         effectID 
sigma^2.2  0.0120  0.1093    203     no          paperID 
sigma^2.3  0.0276  0.1662    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 678) = 389036074.0786, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 678) = 0.6097, p-val = 0.6089

Model Results:

                        estimate      se     tval   df    pval    ci.lb   ci.ub 
dispersal_typeB          -0.0918  0.0756  -1.2150  678  0.2248  -0.2402  0.0566 
dispersal_typebreeding   -0.0125  0.0445  -0.2807  678  0.7790  -0.0998  0.0749 
dispersal_typenatal      -0.0262  0.0330  -0.7939  678  0.4275  -0.0910  0.0386 
                          
dispersal_typeB           
dispersal_typebreeding    
dispersal_typenatal       

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_disp)

   R2_marginal R2_conditional 
   0.002050319    0.266623410 

orchard_plot(mod_disp, mod = "dispersal_type", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S5. Mean fitness effect of natal and breeding dispersal. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (blue circles) scaled by precision (1/SE).

tab_disp <- all_models(mod_disp,  mod = "dispersal_type", type = "normal")

# saving tab_disp as RDS
saveRDS(tab_disp, here("Rdata", "tab_disp.RDS"))

Table S6. Mean estimates of natal and breeding dispersal with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

tab_disp <-readRDS(here("Rdata", "tab_disp.RDS"))

tab_disp

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
B	-0.092	-0.241	0.057	0.226	-0.867	0.683
Breeding	-0.013	-0.100	0.075	0.772	-0.779	0.753
Natal	-0.027	-0.092	0.038	0.410	-0.791	0.736
B-Breeding	0.079	-0.084	0.242	0.343	-0.699	0.857
B-Natal	0.065	-0.087	0.216	0.403	-0.711	0.840
Breeding-Natal	-0.014	-0.100	0.071	0.742	-0.780	0.751

Dispersal phase

mod_phase <- rma.mv(yi = yi, V = VCV,
               mod = ~ dispersal_phase - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_phase)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.3467   556.6933   566.6933   589.2964   566.7825   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1093  0.3306    681     no         effectID 
sigma^2.2  0.0095  0.0973    203     no          paperID 
sigma^2.3  0.0316  0.1777    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 679) = 389782918.7797, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 679) = 1.0351, p-val = 0.3558

Model Results:

                            estimate      se     tval   df    pval    ci.lb 
dispersal_phaseprospection   -0.0664  0.0462  -1.4369  679  0.1512  -0.1572 
dispersal_phasesettlement    -0.0209  0.0318  -0.6558  679  0.5122  -0.0833 
                             ci.ub    
dispersal_phaseprospection  0.0244    
dispersal_phasesettlement   0.0416    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_phase)

   R2_marginal R2_conditional 
   0.002094598    0.274606339 

orchard_plot(mod_phase, mod = "dispersal_phase", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S6. Mean fitness effect of prospection and settlement dispersal phases. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (blue circles) scaled by precision (1/SE).

tab_phase <- all_models(mod_phase,  mod = "dispersal_phase", type = "normal")

# saving tab_phase as RDS
saveRDS(tab_phase, here("Rdata", "tab_phase.RDS"))

Table S7. Mean estimates of prospection and settlement dispersal phases with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

tab_phase <-readRDS(here("Rdata", "tab_phase.RDS"))

tab_phase

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
Prospection	-0.067	-0.158	0.025	0.152	-0.836	0.703
Settlement	-0.022	-0.085	0.041	0.493	-0.788	0.745
Prospection-Settlement	0.045	-0.042	0.132	0.314	-0.724	0.813

Sex

mod_sex <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.1187   556.2373   568.2373   595.3522   568.3625   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1091  0.3304    681     no         effectID 
sigma^2.2  0.0119  0.1091    203     no          paperID 
sigma^2.3  0.0288  0.1696    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 678) = 389669936.4341, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 678) = 0.9185, p-val = 0.4315

Model Results:

      estimate      se     tval   df    pval    ci.lb   ci.ub    
sexB   -0.0734  0.0447  -1.6415  678  0.1012  -0.1613  0.0144    
sexF   -0.0116  0.0352  -0.3305  678  0.7412  -0.0807  0.0574    
sexM   -0.0158  0.0364  -0.4341  678  0.6644  -0.0873  0.0557    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_sex)

   R2_marginal R2_conditional 
   0.003412558    0.274010406 

orchard_plot(mod_sex, mod = "sex", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S7. Mean fitness effect of females, males, or both sexes grouped. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (blue circles) scaled by precision (1/SE).

tab_sex <- all_models(mod_sex,  mod = "sex", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_sex, here("Rdata", "tab_sex.RDS"))

Table S8. Mean estimates of females, males, or both sexes grouped with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

tab_sex <-readRDS(here("Rdata", "tab_sex.RDS"))

tab_sex

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
B	-0.073	-0.162	0.015	0.102	-0.841	0.694
F	-0.013	-0.083	0.056	0.711	-0.779	0.753
M	-0.017	-0.088	0.055	0.650	-0.782	0.749
B-F	0.060	-0.032	0.152	0.199	-0.708	0.828
B-M	0.057	-0.037	0.151	0.237	-0.711	0.825
F-M	-0.003	-0.062	0.055	0.908	-0.768	0.761

Life history stage

3.1 Higher life history grouping

mod_age1 <- rma.mv(yi = yi,
               V = VCV,
               mod = ~ age_class_clean - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_age1)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-276.0583   552.1165   564.1165   591.2314   564.2417   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1094  0.3308    681     no         effectID 
sigma^2.2  0.0087  0.0934    203     no          paperID 
sigma^2.3  0.0304  0.1743    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 678) = 386758443.0763, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 678) = 2.2665, p-val = 0.0796

Model Results:

                         estimate      se     tval   df    pval    ci.lb 
age_class_cleanadult      -0.0268  0.0315  -0.8488  678  0.3963  -0.0887 
age_class_cleanjuvenile    0.0221  0.0506   0.4373  678  0.6620  -0.0772 
age_class_cleanmix        -0.1562  0.0653  -2.3925  678  0.0170  -0.2844 
                           ci.ub    
age_class_cleanadult      0.0352    
age_class_cleanjuvenile   0.1215    
age_class_cleanmix       -0.0280  * 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_age1)

   R2_marginal R2_conditional 
    0.00926677     0.27015532 

orchard_plot(mod_age1, mod = "age_class_clean", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S8. Mean fitness effect of life history stage (adult, juvenile, or mixed). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (blue circles) scaled by precision (1/SE).

tab_age1 <- all_models(mod_age1,  mod = "age_class_clean", type = "normal")

# saving tab_age1 as RDS
saveRDS(tab_age1, here("Rdata", "tab_age1.RDS"))

Table S9. Mean estimates of life history stage (adult, juvenile, or mixed) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

tab_age1 <-readRDS(here("Rdata", "tab_age1.RDS"))

tab_age1

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
Adult	-0.004	-0.061	0.052	0.879	-0.735	0.726
Juvenile	0.012	-0.080	0.105	0.794	-0.722	0.747
Mix	-0.101	-0.224	0.022	0.108	-0.839	0.638
Adult-Juvenile	0.017	-0.072	0.106	0.713	-0.717	0.751
Adult-Mix	-0.096	-0.217	0.025	0.118	-0.835	0.642
Juvenile-Mix	-0.113	-0.254	0.027	0.115	-0.855	0.629

3.2 Finer life stage grouping

# age_class

mod_age2 <- rma.mv(yi = yi,
               V = VCV,
               mod = ~ age_class - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_age2)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-274.5502   549.1003   569.1003   614.2326   569.4321   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1099  0.3315    681     no         effectID 
sigma^2.2  0.0090  0.0950    203     no          paperID 
sigma^2.3  0.0297  0.1724    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 674) = 384503968.8751, p-val < .0001

Test of Moderators (coefficients 1:7):
F(df1 = 7, df2 = 674) = 1.2924, p-val = 0.2512

Model Results:

              estimate      se     tval   df    pval    ci.lb    ci.ub    
age_classA     -0.0194  0.0345  -0.5627  674  0.5738  -0.0870   0.0483    
age_classJ      0.0190  0.0512   0.3708  674  0.7109  -0.0815   0.1194    
age_classJA    -0.2980  0.1829  -1.6288  674  0.1038  -0.6571   0.0612    
age_classJY    -0.2096  0.0973  -2.1542  674  0.0316  -0.4006  -0.0186  * 
age_classJYA   -0.0777  0.0923  -0.8414  674  0.4004  -0.2590   0.1036    
age_classY     -0.0451  0.0610  -0.7387  674  0.4604  -0.1648   0.0747    
age_classYA    -0.0523  0.0495  -1.0570  674  0.2909  -0.1495   0.0449    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_age2)

   R2_marginal R2_conditional 
    0.01393744     0.27099634 

orchard_plot(mod_age2, mod = "age_class", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S9. Mean fitness effect of life history stage at the level of precision reported (A: adult, Y: yearling, J: juvenile). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (blue circles) scaled by precision (1/SE).

# table not created for this (not meaningful for some levels)

Fitness metric

3.3 Higher level

mod_fit1 <- rma.mv(yi = yi, 
               V = VCV,
               mod = ~ fitness_higher_level - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_fit1)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.1028   556.2057   568.2057   595.3206   568.3309   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1095  0.3309    681     no         effectID 
sigma^2.2  0.0112  0.1060    203     no          paperID 
sigma^2.3  0.0285  0.1689    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 678) = 389732718.8287, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 678) = 0.9513, p-val = 0.4153

Model Results:

                                      estimate      se     tval   df    pval 
fitness_higher_levelindirect fitness   -0.2364  0.1700  -1.3908  678  0.1647 
fitness_higher_levelreproduction       -0.0177  0.0333  -0.5326  678  0.5945 
fitness_higher_levelsurvival           -0.0369  0.0344  -1.0734  678  0.2835 
                                        ci.lb   ci.ub    
fitness_higher_levelindirect fitness  -0.5702  0.0973    
fitness_higher_levelreproduction      -0.0831  0.0476    
fitness_higher_levelsurvival          -0.1043  0.0306    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_fit1)

   R2_marginal R2_conditional 
   0.002751919    0.268487953 

orchard_plot(mod_fit1, mod = "fitness_higher_level", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S10. Mean fitness effect of different fitness components (survival, reproduction, and indirect fitness). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (blue circles) scaled by precision (1/SE).

tab_fit1 <- all_models(mod_fit1,  mod = "fitness_higher_level", type = "normal")

# saving tab_fit1 as RDS
saveRDS(tab_fit1, here("Rdata", "tab_fit1.RDS"))

Table S10. Mean estimates of different fitness components (survival, reproduction, and indirect fitness) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

tab_fit1 <-readRDS(here("Rdata", "tab_fit1.RDS"))

tab_fit1

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
Indirect fitness	-0.237	-0.572	0.098	0.166	-1.068	0.595
Reproduction	-0.019	-0.084	0.047	0.580	-0.782	0.745
Survival	-0.038	-0.106	0.030	0.273	-0.802	0.726
Indirect fitness-Reproduction	0.218	-0.117	0.554	0.202	-0.613	1.050
Indirect fitness-Survival	0.199	-0.136	0.534	0.244	-0.632	1.030
Reproduction-Survival	-0.019	-0.079	0.040	0.525	-0.783	0.744

3.4 Finer level

• need to think what to do with classes less than 5 studies

mod_fit2 <- rma.mv(yi = yi, 
               V = VCV,
               mod = ~ fitness_metric_clean - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_fit2)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-277.1187   554.2374   576.2374   625.8666   576.6368   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1096  0.3310    681     no         effectID 
sigma^2.2  0.0086  0.0930    203     no          paperID 
sigma^2.3  0.0336  0.1833    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 673) = 387473739.5403, p-val < .0001

Test of Moderators (coefficients 1:8):
F(df1 = 8, df2 = 673) = 0.6888, p-val = 0.7016

Model Results:

                                                           estimate      se 
fitness_metric_cleanage at first reproduction               -0.0349  0.1056 
fitness_metric_cleanindirect fitness                        -0.2309  0.1709 
fitness_metric_cleanlifetime breeding success               -0.0056  0.0388 
fitness_metric_cleanlifetime reproductive success           -0.0322  0.0425 
fitness_metric_cleanoffspring reproduction                   0.0387  0.1065 
fitness_metric_cleanoffspring survival                       0.0029  0.0480 
fitness_metric_cleanreproductive lifespan and/or attempts   -0.0671  0.1021 
fitness_metric_cleansurvival                                -0.0572  0.0382 
                                                              tval   df    pval 
fitness_metric_cleanage at first reproduction              -0.3306  673  0.7411 
fitness_metric_cleanindirect fitness                       -1.3510  673  0.1772 
fitness_metric_cleanlifetime breeding success              -0.1450  673  0.8847 
fitness_metric_cleanlifetime reproductive success          -0.7588  673  0.4482 
fitness_metric_cleanoffspring reproduction                  0.3629  673  0.7168 
fitness_metric_cleanoffspring survival                      0.0612  673  0.9512 
fitness_metric_cleanreproductive lifespan and/or attempts  -0.6573  673  0.5112 
fitness_metric_cleansurvival                               -1.4968  673  0.1349 
                                                             ci.lb   ci.ub    
fitness_metric_cleanage at first reproduction              -0.2422  0.1724    
fitness_metric_cleanindirect fitness                       -0.5665  0.1047    
fitness_metric_cleanlifetime breeding success              -0.0817  0.0705    
fitness_metric_cleanlifetime reproductive success          -0.1157  0.0512    
fitness_metric_cleanoffspring reproduction                 -0.1705  0.2478    
fitness_metric_cleanoffspring survival                     -0.0912  0.0971    
fitness_metric_cleanreproductive lifespan and/or attempts  -0.2676  0.1334    
fitness_metric_cleansurvival                               -0.1323  0.0178    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_fit2)

   R2_marginal R2_conditional 
   0.006389204    0.282906100 

orchard_plot(mod_fit2, mod = "fitness_metric_clean", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S11. Mean fitness effect of different fitness components at a finer resolution. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (blue circles) scaled by precision (1/SE).

# table not created for this (not meaningful for some levels)

Multi-generation

mod_gen <- rma.mv(yi = yi, 
                V = VCV,
               mod = ~ whose_fitness - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_gen)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.4694   556.9388   566.9388   589.5419   567.0280   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1089  0.3300    681     no         effectID 
sigma^2.2  0.0087  0.0934    203     no          paperID 
sigma^2.3  0.0339  0.1842    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 679) = 390484843.0291, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 679) = 1.1121, p-val = 0.3295

Model Results:

                         estimate      se     tval   df    pval    ci.lb 
whose_fitnessdescendant    0.0058  0.0462   0.1259  679  0.8999  -0.0849 
whose_fitnessfocal        -0.0374  0.0312  -1.1998  679  0.2306  -0.0986 
                          ci.ub    
whose_fitnessdescendant  0.0966    
whose_fitnessfocal       0.0238    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_gen)

   R2_marginal R2_conditional 
   0.001680206    0.282516227 

orchard_plot(mod_gen, mod = "whose_fitness", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S12. Mean fitness effect of non-dispersers and dispersers and their descendants. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (blue circles) scaled by precision (1/SE).

tab_gen <- all_models(mod_gen,  mod = "whose_fitness", type = "normal")

# saving tab_gen as RDS
saveRDS(tab_gen, here("Rdata", "tab_gen.RDS"))

Table S11. Mean estimates of non-dispersers and dispersers and their descendants with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

tab_gen <-readRDS(here("Rdata", "tab_gen.RDS"))

tab_gen

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
Individual	-0.038	-0.100	0.023	0.219	-0.808	0.731
Offspring	0.005	-0.086	0.096	0.907	-0.767	0.777
Individual-Offspring	0.044	-0.035	0.123	0.275	-0.727	0.815

Study duration

3.5 Normal analysis

# getting absolute values for effect sizes as we expect shorter years will have more fluctuations
dat$ln_study_duration <- log(dat$study_duration+1)

mod_dur <- rma.mv(yi = yi, 
                V = VCV,
               mod = ~ ln_study_duration,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_dur)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.7049   557.4098   567.4098   590.0129   567.4989   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1095  0.3308    681     no         effectID 
sigma^2.2  0.0113  0.1064    203     no          paperID 
sigma^2.3  0.0291  0.1705    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 679) = 390307264.2821, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 679) = 0.0244, p-val = 0.8759

Model Results:

                   estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt             -0.0215  0.0593  -0.3622  679  0.7173  -0.1378  0.0949    
ln_study_duration   -0.0034  0.0215  -0.1562  679  0.8759  -0.0455  0.0388    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_dur)

   R2_marginal R2_conditional 
  0.0000656924   0.2696351007 

bubble_plot(mod_dur,
                mod = "ln_study_duration",
                group = "paperID",
                xlab = "Study duration",
                ylab = "Effect Size: Zr",
                       g = TRUE)

Figure S13. Add text.

# absolute value analysis
mod_dur_abs <- rma.mv(yi = yi_abs,
                V = VCVa,
               mod = ~ ln_study_duration,
               data = dat,
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_dur_abs)

Multivariate Meta-Analysis Model (k = 681; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-12.4516   24.9032   34.9032   57.5063   34.9923   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0481  0.2193    681     no         effectID 
sigma^2.2  0.0073  0.0853    203     no          paperID 
sigma^2.3  0.0128  0.1130    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 679) = 334226677.4018, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 679) = 2.7486, p-val = 0.0978

Model Results:

                   estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt             -0.0149  0.0411  -0.3625  679  0.7171  -0.0956  0.0658    
ln_study_duration    0.0246  0.0148   1.6579  679  0.0978  -0.0045  0.0536  . 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_dur_abs)

   R2_marginal R2_conditional 
   0.007677259    0.299605809 

bubble_plot(mod_dur_abs,
                mod = "ln_study_duration",
                group = "paperID",
                xlab = "Study duration",
                ylab = "Effect Size: Zr",
                       g = TRUE)

Figure S14. Add text.

Selection bias

3.6 Normal analysis

mod_focus <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ fitness_main_focus - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                sparse = TRUE)
summary(mod_focus)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.7221   557.4441   567.4441   590.0472   567.5333   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1093  0.3306    681     no         effectID 
sigma^2.2  0.0116  0.1079    203     no          paperID 
sigma^2.3  0.0291  0.1705    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 679) = 390349329.8724, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 679) = 0.5257, p-val = 0.5914

Model Results:

                     estimate      se     tval   df    pval    ci.lb   ci.ub    
fitness_main_focusN   -0.0181  0.0463  -0.3913  679  0.6957  -0.1090  0.0728    
fitness_main_focusY   -0.0326  0.0318  -1.0253  679  0.3056  -0.0951  0.0298    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_focus)

   R2_marginal R2_conditional 
  0.0002547069   0.2715288977 

orchard_plot(mod_focus, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S17. Add text.

tab_focus <- all_models(mod_focus,  mod = "fitness_main_focus", type = "normal")

# saving tab_focus as RDS
saveRDS(tab_focus, here("Rdata", "tab_focus.RDS"))

Table S14. Add text.

tab_focus <-readRDS(here("Rdata", "tab_focus.RDS"))

tab_focus

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
N	-0.019	-0.110	0.073	0.690	-0.787	0.750
Y	-0.034	-0.096	0.029	0.292	-0.799	0.732
N-Y	-0.015	-0.103	0.073	0.734	-0.783	0.752

3.7 Comparing normal model to heteroscedastic model

# homo
mod_focus1 <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ fitness_main_focus - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                method = "ML",
                sparse = TRUE)
summary(mod_focus1)

Multivariate Meta-Analysis Model (k = 681; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-279.6163  9646.0635   569.2326   591.8504   569.3215   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1095  0.3309    681     no         effectID 
sigma^2.2  0.0114  0.1068    203     no          paperID 
sigma^2.3  0.0278  0.1668    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 679) = 390349329.8724, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 679) = 0.5078, p-val = 0.6020

Model Results:

                     estimate      se     tval   df    pval    ci.lb   ci.ub    
fitness_main_focusN   -0.0179  0.0460  -0.3890  679  0.6974  -0.1082  0.0724    
fitness_main_focusY   -0.0318  0.0315  -1.0078  679  0.3139  -0.0937  0.0301    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_focus1)

   R2_marginal R2_conditional 
  0.0002358119   0.2638866485 

# hetero
mod_focus2 <- rma.mv(yi = yi, 
                 V = VCV,
                 mod = ~ fitness_main_focus,
                 data = dat, 
                 random = list(
                   ~ fitness_main_focus | effectID,
                   ~ 1 | paperID,
                   ~ 1 | species_cleaned),
                 struct = "DIAG",
                 method = "ML",
                 test = "t",
                 sparse = TRUE,
                  control=list(optimizer="nlminb", rel.tol=1e-8))
summary(mod_focus2)

Multivariate Meta-Analysis Model (k = 681; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-279.4838  9645.7985   570.9676   598.1089   571.0922   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0116  0.1075    203     no          paperID 
sigma^2.2  0.0271  0.1646    147     no  species_cleaned 

outer factor: effectID           (nlvls = 681)
inner factor: fitness_main_focus (nlvls = 2)

            estim    sqrt  k.lvl  fixed  level 
tau^2.1    0.1035  0.3217    162     no      N 
tau^2.2    0.1115  0.3340    519     no      Y 

Test for Residual Heterogeneity:
QE(df = 679) = 390349329.8724, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 679) = 0.0995, p-val = 0.7525

Model Results:

                     estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt               -0.0173  0.0453  -0.3822  679  0.7024  -0.1063  0.0717    
fitness_main_focusY   -0.0138  0.0438  -0.3155  679  0.7525  -0.0999  0.0722    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#r2_ml(mod_focusA1)

# they are not significantly different
anova(mod_focus1, mod_focus2)

        df      AIC      BIC     AICc    logLik    LRT   pval             QE 
Full     6 570.9676 598.1089 571.0922 -279.4838               390349329.8724 
Reduced  5 569.2326 591.8504 569.3215 -279.6163 0.2650 0.6067 390349329.8724 

# homo
orchard_plot(mod_focus1, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S18. Add text.

# hetero
orchard_plot(mod_focus2, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S19. Add text.

Confirmation bias

3.8 Normal analysis

mod_confirm <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ confirmation_bias - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                sparse = TRUE)
summary(mod_confirm)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.1402   556.2805   570.2805   601.9042   570.4479   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1095  0.3309    681     no         effectID 
sigma^2.2  0.0170  0.1306    203     no          paperID 
sigma^2.3  0.0220  0.1483    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 677) = 390198859.2817, p-val < .0001

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 677) = 0.3307, p-val = 0.8574

Model Results:

                          estimate      se     tval   df    pval    ci.lb 
confirmation_biasBenefit    0.0109  0.0842   0.1289  677  0.8975  -0.1545 
confirmation_biasCost      -0.0395  0.0402  -0.9829  677  0.3260  -0.1183 
confirmation_biasMixed     -0.0290  0.0401  -0.7235  677  0.4696  -0.1077 
confirmation_biasNeutral   -0.0099  0.0473  -0.2086  677  0.8348  -0.1028 
                           ci.ub    
confirmation_biasBenefit  0.1763    
confirmation_biasCost     0.0394    
confirmation_biasMixed    0.0497    
confirmation_biasNeutral  0.0831    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_confirm)

   R2_marginal R2_conditional 
   0.001217654    0.263662784 

orchard_plot(mod_confirm, mod = "confirmation_bias", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S20. Add text.

tab_confirm <- all_models(mod_confirm,  mod = "confirmation_bias", type = "normal")

# saving tab_confirm as RDS
saveRDS(tab_confirm, here("Rdata", "tab_confirm.RDS"))

Table S15. Add text.

tab_confirm <-readRDS(here("Rdata", "tab_confirm.RDS"))

tab_confirm

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
Benefit	0.012	-0.154	0.178	0.884	-0.764	0.788
Cost	-0.042	-0.122	0.037	0.296	-0.804	0.720
Mixed	-0.028	-0.107	0.050	0.483	-0.790	0.734
Neutral	-0.009	-0.102	0.084	0.852	-0.772	0.755
Benefit-Cost	-0.055	-0.227	0.118	0.536	-0.832	0.723
Benefit-Mixed	-0.040	-0.212	0.132	0.645	-0.818	0.737
Benefit-Neutral	-0.021	-0.198	0.156	0.815	-0.800	0.757
Cost-Mixed	0.014	-0.081	0.110	0.770	-0.750	0.778
Cost-Neutral	0.033	-0.071	0.138	0.531	-0.732	0.799
Mixed-Neutral	0.019	-0.085	0.124	0.719	-0.746	0.784

3.9 Comparing normal model to heteroscedastic model

# homo
mod_confirm1 <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ confirmation_bias - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                method = "ML",
                sparse = TRUE)
summary(mod_confirm1)

Multivariate Meta-Analysis Model (k = 681; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-279.4495  9645.7299   572.8989   604.5638   573.0653   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1099  0.3314    681     no         effectID 
sigma^2.2  0.0156  0.1250    203     no          paperID 
sigma^2.3  0.0209  0.1446    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 677) = 390198859.2817, p-val < .0001

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 677) = 0.3192, p-val = 0.8652

Model Results:

                          estimate      se     tval   df    pval    ci.lb 
confirmation_biasBenefit    0.0105  0.0832   0.1262  677  0.8996  -0.1528 
confirmation_biasCost      -0.0385  0.0396  -0.9719  677  0.3315  -0.1163 
confirmation_biasMixed     -0.0276  0.0394  -0.7004  677  0.4839  -0.1050 
confirmation_biasNeutral   -0.0093  0.0466  -0.1999  677  0.8416  -0.1009 
                           ci.ub    
confirmation_biasBenefit  0.1738    
confirmation_biasCost     0.0393    
confirmation_biasMixed    0.0498    
confirmation_biasNeutral  0.0822    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_confirm1)

   R2_marginal R2_conditional 
   0.001178183    0.250425872 

# hetero
mod_confirm2 <- rma.mv(yi = yi, 
                 V = VCV,
                 mod = ~ confirmation_bias,
                 data = dat, 
                 random = list(
                   ~ confirmation_bias | effectID,
                   ~ 1 | paperID,
                   ~ 1 | species_cleaned),
                 struct = "DIAG",
                 method = "ML",
                 test = "t",
                 sparse = TRUE,
                  control=list(optimizer="nlminb", rel.tol=1e-8))
summary(mod_confirm2)

Multivariate Meta-Analysis Model (k = 681; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-241.1711  9569.1732   502.3422   547.5778   502.6706   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0149  0.1221    203     no          paperID 
sigma^2.2  0.0133  0.1155    147     no  species_cleaned 

outer factor: effectID          (nlvls = 681)
inner factor: confirmation_bias (nlvls = 4)

            estim    sqrt  k.lvl  fixed    level 
tau^2.1    0.0295  0.1718     32     no  Benefit 
tau^2.2    0.1324  0.3638    210     no     Cost 
tau^2.3    0.0633  0.2516    289     no    Mixed 
tau^2.4    0.2045  0.4523    150     no  Neutral 

Test for Residual Heterogeneity:
QE(df = 677) = 390198859.2817, p-val < .0001

Test of Moderators (coefficients 2:4):
F(df1 = 3, df2 = 677) = 0.4440, p-val = 0.7216

Model Results:

                          estimate      se     tval   df    pval    ci.lb 
intrcpt                     0.0230  0.0594   0.3869  677  0.6990  -0.0937 
confirmation_biasCost      -0.0730  0.0648  -1.1278  677  0.2598  -0.2002 
confirmation_biasMixed     -0.0595  0.0617  -0.9637  677  0.3355  -0.1807 
confirmation_biasNeutral   -0.0467  0.0716  -0.6525  677  0.5143  -0.1873 
                           ci.ub    
intrcpt                   0.1396    
confirmation_biasCost     0.0541    
confirmation_biasMixed    0.0617    
confirmation_biasNeutral  0.0939    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# they are not significantly different
anova(mod_confirm1, mod_confirm2)

        df      AIC      BIC     AICc    logLik     LRT   pval             QE 
Full    10 502.3422 547.5778 502.6706 -241.1711                390198859.2817 
Reduced  7 572.8989 604.5638 573.0653 -279.4495 76.5567 <.0001 390198859.2817 

# homo
orchard_plot(mod_confirm1, mod = "confirmation_bias", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S21. Add text.

# hetero
orchard_plot(mod_confirm2, mod = "confirmation_bias", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S22. Add text.

4 Interaction meta-regression models (2 moderators)

Sex x Species class

# the interaction between 2 categorical variables are the same as creating a new variable with 2 categories combined so we use that approach 

dat$sex_species_class <- as.factor(paste(dat$sex, dat$species_class ,sep = "_"))

mod_sex_species_class <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_species_class - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex_species_class)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-272.8754   545.7507   581.7507   662.7739   582.8079   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1093  0.3306    681     no         effectID 
sigma^2.2  0.0113  0.1065    203     no          paperID 
sigma^2.3  0.0309  0.1757    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 666) = 383462543.2723, p-val < .0001

Test of Moderators (coefficients 1:15):
F(df1 = 15, df2 = 666) = 0.7720, p-val = 0.7096

Model Results:

                                   estimate      se     tval   df    pval 
sex_species_classB_Aves             -0.0009  0.0582  -0.0161  666  0.9872 
sex_species_classB_Insecta          -0.1461  0.2753  -0.5308  666  0.5958 
sex_species_classB_Mammalia         -0.1785  0.0671  -2.6597  666  0.0080 
sex_species_classB_Reptilia          0.1585  0.2731   0.5803  666  0.5619 
sex_species_classF_Actinopterygii   -0.0682  0.2069  -0.3295  666  0.7419 
sex_species_classF_Arachnida         0.0406  0.2810   0.1446  666  0.8851 
sex_species_classF_Aves             -0.0360  0.0443  -0.8122  666  0.4170 
sex_species_classF_Insecta           0.2396  0.1865   1.2847  666  0.1994 
sex_species_classF_Mammalia         -0.0140  0.0528  -0.2649  666  0.7912 
sex_species_classF_Reptilia          0.0719  0.1565   0.4594  666  0.6461 
sex_species_classM_Actinopterygii    0.0688  0.2063   0.3335  666  0.7388 
sex_species_classM_Aves             -0.0345  0.0447  -0.7729  666  0.4398 
sex_species_classM_Insecta          -0.1425  0.3778  -0.3771  666  0.7062 
sex_species_classM_Mammalia         -0.0134  0.0571  -0.2354  666  0.8140 
sex_species_classM_Reptilia         -0.0233  0.2318  -0.1005  666  0.9200 
                                     ci.lb    ci.ub     
sex_species_classB_Aves            -0.1152   0.1134     
sex_species_classB_Insecta         -0.6867   0.3945     
sex_species_classB_Mammalia        -0.3104  -0.0467  ** 
sex_species_classB_Reptilia        -0.3777   0.6947     
sex_species_classF_Actinopterygii  -0.4743   0.3380     
sex_species_classF_Arachnida       -0.5111   0.5923     
sex_species_classF_Aves            -0.1230   0.0510     
sex_species_classF_Insecta         -0.1266   0.6058     
sex_species_classF_Mammalia        -0.1176   0.0896     
sex_species_classF_Reptilia        -0.2354   0.3792     
sex_species_classM_Actinopterygii  -0.3363   0.4739     
sex_species_classM_Aves            -0.1222   0.0532     
sex_species_classM_Insecta         -0.8842   0.5993     
sex_species_classM_Mammalia        -0.1256   0.0987     
sex_species_classM_Reptilia        -0.4784   0.4318     

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_sex_species_class)

   R2_marginal R2_conditional 
    0.01839299     0.29193461 

orchard_plot(mod_sex_species_class, mod = "sex_species_class", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S23. Add text.

tab_sex_species_class <- all_models(mod_sex_species_class,  mod = "sex_species_class", type = "normal")

# saving tab_sex_species_class as RDS
saveRDS(tab_sex_species_class, here("Rdata", "tab_sex_species_class.RDS"))

Table S16. Add text.

tab_sex_species_class <-readRDS(here("Rdata", "tab_sex_species_class.RDS"))

tab_sex_species_class

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
B_Aves	-0.001	-0.116	0.114	0.987	-0.777	0.775
B_Insecta	-0.146	-0.689	0.397	0.597	-1.086	0.794
B_Mammalia	-0.180	-0.312	-0.047	0.008	-0.959	0.599
B_Reptilia	0.158	-0.381	0.696	0.565	-0.780	1.095
F_Actinopterygii	-0.069	-0.477	0.339	0.741	-0.938	0.800
F_Arachnida	0.040	-0.514	0.594	0.887	-0.906	0.987
F_Aves	-0.036	-0.124	0.051	0.417	-0.809	0.736
F_Insecta	0.240	-0.128	0.607	0.201	-0.611	1.091
F_Mammalia	-0.019	-0.124	0.087	0.731	-0.793	0.756
F_Reptilia	0.071	-0.238	0.380	0.651	-0.756	0.898
M_Actinopterygii	0.068	-0.339	0.475	0.742	-0.800	0.937
M_Aves	-0.035	-0.123	0.054	0.441	-0.807	0.738
M_Insecta	-0.143	-0.888	0.602	0.707	-1.212	0.927
M_Mammalia	-0.016	-0.129	0.097	0.778	-0.792	0.760
M_Reptilia	-0.024	-0.481	0.433	0.918	-0.917	0.869
B_Aves-B_Insecta	-0.145	-0.700	0.409	0.607	-1.092	0.801
B_Aves-B_Mammalia	-0.179	-0.346	-0.012	0.036	-0.964	0.607
B_Aves-B_Reptilia	0.159	-0.388	0.706	0.569	-0.784	1.101
B_Aves-F_Actinopterygii	-0.068	-0.489	0.353	0.752	-0.943	0.807
B_Aves-F_Arachnida	0.041	-0.523	0.605	0.886	-0.912	0.994
B_Aves-F_Aves	-0.035	-0.165	0.094	0.594	-0.814	0.743
B_Aves-F_Insecta	0.241	-0.141	0.623	0.217	-0.617	1.098
B_Aves-F_Mammalia	-0.018	-0.165	0.130	0.815	-0.799	0.764
B_Aves-F_Reptilia	0.072	-0.253	0.396	0.663	-0.761	0.905
B_Aves-M_Actinopterygii	0.069	-0.351	0.490	0.746	-0.806	0.944
B_Aves-M_Aves	-0.034	-0.163	0.096	0.611	-0.812	0.745
B_Aves-M_Insecta	-0.142	-0.894	0.611	0.712	-1.217	0.933
B_Aves-M_Mammalia	-0.015	-0.169	0.138	0.845	-0.798	0.767
B_Aves-M_Reptilia	-0.023	-0.491	0.445	0.923	-0.922	0.876
B_Insecta-B_Mammalia	-0.033	-0.591	0.524	0.906	-0.982	0.915
B_Insecta-B_Reptilia	0.304	-0.459	1.067	0.435	-0.779	1.386
B_Insecta-F_Actinopterygii	0.077	-0.601	0.756	0.823	-0.947	1.102
B_Insecta-F_Arachnida	0.186	-0.589	0.962	0.637	-0.904	1.277
B_Insecta-F_Aves	0.110	-0.439	0.659	0.694	-0.834	1.054
B_Insecta-F_Insecta	0.386	-0.269	1.041	0.248	-0.623	1.395
B_Insecta-F_Mammalia	0.128	-0.424	0.680	0.650	-0.818	1.073
B_Insecta-F_Reptilia	0.217	-0.406	0.841	0.494	-0.771	1.206
B_Insecta-M_Actinopterygii	0.215	-0.463	0.892	0.534	-0.809	1.239
B_Insecta-M_Aves	0.112	-0.438	0.661	0.690	-0.832	1.056
B_Insecta-M_Insecta	0.004	-0.918	0.925	0.994	-1.196	1.203
B_Insecta-M_Mammalia	0.130	-0.423	0.684	0.645	-0.816	1.076
B_Insecta-M_Reptilia	0.122	-0.587	0.831	0.735	-0.922	1.167
B_Mammalia-B_Reptilia	0.337	-0.210	0.885	0.227	-0.605	1.280
B_Mammalia-F_Actinopterygii	0.111	-0.312	0.534	0.607	-0.765	0.987
B_Mammalia-F_Arachnida	0.220	-0.347	0.786	0.446	-0.734	1.174
B_Mammalia-F_Aves	0.144	-0.004	0.291	0.056	-0.638	0.925
B_Mammalia-F_Insecta	0.419	0.035	0.804	0.033	-0.439	1.278
B_Mammalia-F_Mammalia	0.161	0.020	0.303	0.026	-0.619	0.942
B_Mammalia-F_Reptilia	0.251	-0.075	0.577	0.132	-0.583	1.085
B_Mammalia-M_Actinopterygii	0.248	-0.175	0.671	0.250	-0.628	1.124
B_Mammalia-M_Aves	0.145	-0.002	0.293	0.054	-0.636	0.927
B_Mammalia-M_Insecta	0.037	-0.717	0.791	0.923	-1.039	1.113
B_Mammalia-M_Mammalia	0.163	0.015	0.312	0.031	-0.618	0.945
B_Mammalia-M_Reptilia	0.156	-0.314	0.625	0.515	-0.744	1.055
B_Reptilia-F_Actinopterygii	-0.226	-0.897	0.444	0.507	-1.245	0.793
B_Reptilia-F_Arachnida	-0.117	-0.887	0.652	0.764	-1.204	0.969
B_Reptilia-F_Aves	-0.194	-0.735	0.347	0.482	-1.133	0.745
B_Reptilia-F_Insecta	0.082	-0.565	0.729	0.803	-0.922	1.086
B_Reptilia-F_Mammalia	-0.176	-0.718	0.366	0.524	-1.116	0.764
B_Reptilia-F_Reptilia	-0.087	-0.653	0.479	0.764	-1.040	0.867
B_Reptilia-M_Actinopterygii	-0.089	-0.760	0.581	0.794	-1.109	0.930
B_Reptilia-M_Aves	-0.192	-0.733	0.349	0.486	-1.131	0.747
B_Reptilia-M_Insecta	-0.300	-1.217	0.616	0.520	-1.496	0.895
B_Reptilia-M_Mammalia	-0.174	-0.718	0.370	0.530	-1.115	0.767
B_Reptilia-M_Reptilia	-0.182	-0.860	0.497	0.599	-1.206	0.843
F_Actinopterygii-F_Arachnida	0.109	-0.577	0.795	0.755	-0.920	1.138
F_Actinopterygii-F_Aves	0.033	-0.381	0.446	0.877	-0.839	0.904
F_Actinopterygii-F_Insecta	0.308	-0.237	0.854	0.267	-0.633	1.250
F_Actinopterygii-F_Mammalia	0.050	-0.366	0.466	0.812	-0.823	0.923
F_Actinopterygii-F_Reptilia	0.140	-0.366	0.646	0.588	-0.779	1.059
F_Actinopterygii-M_Actinopterygii	0.137	-0.324	0.598	0.559	-0.758	1.032
F_Actinopterygii-M_Aves	0.034	-0.379	0.448	0.871	-0.838	0.906
F_Actinopterygii-M_Insecta	-0.074	-0.921	0.773	0.864	-1.217	1.069
F_Actinopterygii-M_Mammalia	0.053	-0.366	0.471	0.805	-0.821	0.927
F_Actinopterygii-M_Reptilia	0.045	-0.564	0.653	0.885	-0.935	1.024
F_Arachnida-F_Aves	-0.076	-0.635	0.483	0.789	-1.026	0.873
F_Arachnida-F_Insecta	0.200	-0.463	0.862	0.555	-0.814	1.213
F_Arachnida-F_Mammalia	-0.059	-0.620	0.502	0.837	-1.009	0.892
F_Arachnida-F_Reptilia	0.031	-0.600	0.662	0.924	-0.963	1.024
F_Arachnida-M_Actinopterygii	0.028	-0.657	0.714	0.936	-1.001	1.057
F_Arachnida-M_Aves	-0.075	-0.634	0.484	0.793	-1.024	0.875
F_Arachnida-M_Insecta	-0.183	-1.110	0.744	0.699	-1.386	1.021
F_Arachnida-M_Mammalia	-0.056	-0.619	0.506	0.844	-1.008	0.895
F_Arachnida-M_Reptilia	-0.064	-0.780	0.651	0.860	-1.114	0.985
F_Aves-F_Insecta	0.276	-0.098	0.650	0.148	-0.578	1.130
F_Aves-F_Mammalia	0.018	-0.108	0.143	0.782	-0.760	0.795
F_Aves-F_Reptilia	0.107	-0.208	0.422	0.504	-0.722	0.937
F_Aves-M_Actinopterygii	0.104	-0.309	0.518	0.620	-0.767	0.976
F_Aves-M_Aves	0.002	-0.074	0.077	0.967	-0.770	0.773
F_Aves-M_Insecta	-0.107	-0.855	0.642	0.780	-1.179	0.965
F_Aves-M_Mammalia	0.020	-0.112	0.152	0.768	-0.759	0.799
F_Aves-M_Reptilia	0.012	-0.449	0.474	0.959	-0.884	0.908
F_Insecta-F_Mammalia	-0.258	-0.635	0.118	0.179	-1.113	0.597
F_Insecta-F_Reptilia	-0.169	-0.643	0.305	0.485	-1.071	0.733
F_Insecta-M_Actinopterygii	-0.171	-0.717	0.374	0.537	-1.113	0.770
F_Insecta-M_Aves	-0.274	-0.648	0.100	0.151	-1.128	0.580
F_Insecta-M_Insecta	-0.382	-1.165	0.400	0.338	-1.479	0.714
F_Insecta-M_Mammalia	-0.256	-0.635	0.123	0.185	-1.112	0.600
F_Insecta-M_Reptilia	-0.264	-0.846	0.318	0.374	-1.227	0.700
F_Mammalia-F_Reptilia	0.090	-0.228	0.407	0.580	-0.741	0.920
F_Mammalia-M_Actinopterygii	0.087	-0.329	0.503	0.682	-0.786	0.960
F_Mammalia-M_Aves	-0.016	-0.141	0.109	0.802	-0.794	0.762
F_Mammalia-M_Insecta	-0.124	-0.874	0.626	0.745	-1.197	0.949
F_Mammalia-M_Mammalia	0.002	-0.099	0.103	0.965	-0.772	0.776
F_Mammalia-M_Reptilia	-0.006	-0.469	0.458	0.981	-0.902	0.891
F_Reptilia-M_Actinopterygii	-0.003	-0.509	0.503	0.992	-0.922	0.917
F_Reptilia-M_Aves	-0.106	-0.420	0.209	0.511	-0.935	0.724
F_Reptilia-M_Insecta	-0.214	-1.017	0.590	0.602	-1.325	0.898
F_Reptilia-M_Mammalia	-0.087	-0.407	0.233	0.593	-0.919	0.744
F_Reptilia-M_Reptilia	-0.095	-0.542	0.352	0.676	-0.983	0.793
M_Actinopterygii-M_Aves	-0.103	-0.516	0.311	0.625	-0.975	0.769
M_Actinopterygii-M_Insecta	-0.211	-1.058	0.636	0.625	-1.354	0.932
M_Actinopterygii-M_Mammalia	-0.085	-0.503	0.334	0.691	-0.959	0.789
M_Actinopterygii-M_Reptilia	-0.092	-0.701	0.516	0.766	-1.072	0.887
M_Aves-M_Insecta	-0.108	-0.857	0.640	0.777	-1.180	0.964
M_Aves-M_Mammalia	0.018	-0.114	0.151	0.786	-0.761	0.797
M_Aves-M_Reptilia	0.010	-0.451	0.472	0.965	-0.885	0.906
M_Insecta-M_Mammalia	0.126	-0.625	0.878	0.741	-0.948	1.201
M_Insecta-M_Reptilia	0.119	-0.753	0.990	0.789	-1.043	1.280
M_Mammalia-M_Reptilia	-0.008	-0.473	0.457	0.974	-0.905	0.890

Sex x Whose fitness

dat$sex_whose_fitness <- as.factor(paste(dat$sex, dat$whose_fitness ,sep = "_"))

mod_sex_whose_fitness <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_whose_fitness - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex_whose_fitness)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-276.9443   553.8885   571.8885   612.5210   572.1592   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1092  0.3305    681     no         effectID 
sigma^2.2  0.0089  0.0944    203     no          paperID 
sigma^2.3  0.0329  0.1815    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 675) = 387467254.1390, p-val < .0001

Test of Moderators (coefficients 1:6):
F(df1 = 6, df2 = 675) = 0.9701, p-val = 0.4445

Model Results:

                               estimate      se     tval   df    pval    ci.lb 
sex_whose_fitnessB_descendant    0.0204  0.1500   0.1358  675  0.8920  -0.2742 
sex_whose_fitnessB_focal        -0.0783  0.0455  -1.7229  675  0.0854  -0.1676 
sex_whose_fitnessF_descendant   -0.0204  0.0548  -0.3716  675  0.7103  -0.1279 
sex_whose_fitnessF_focal        -0.0096  0.0369  -0.2611  675  0.7941  -0.0820 
sex_whose_fitnessM_descendant    0.0616  0.0629   0.9795  675  0.3277  -0.0619 
sex_whose_fitnessM_focal        -0.0329  0.0380  -0.8651  675  0.3873  -0.1075 
                                ci.ub    
sex_whose_fitnessB_descendant  0.3150    
sex_whose_fitnessB_focal       0.0109  . 
sex_whose_fitnessF_descendant  0.0872    
sex_whose_fitnessF_focal       0.0627    
sex_whose_fitnessM_descendant  0.1851    
sex_whose_fitnessM_focal       0.0417    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_sex_whose_fitness)

   R2_marginal R2_conditional 
   0.006900915    0.282083810 

orchard_plot(mod_sex_whose_fitness, mod = "sex_whose_fitness", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S23. Add text.

tab_sex_whose_fitness <- all_models(mod_sex_whose_fitness,  mod = "sex_whose_fitness", type = "normal")

# saving tab_sex_whose_fitness as RDS
saveRDS(tab_sex_whose_fitness, here("Rdata", "tab_sex_whose_fitness.RDS"))

Table S17. Add text.

tab_sex_whose_fitness <-readRDS(here("Rdata", "tab_sex_whose_fitness.RDS"))

tab_sex_whose_fitness

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
B_individual	-0.078	-0.168	0.011	0.086	-0.849	0.693
B_offspring	0.020	-0.276	0.316	0.894	-0.801	0.841
F_individual	-0.011	-0.084	0.062	0.760	-0.781	0.758
F_offspring	-0.021	-0.129	0.087	0.704	-0.794	0.753
M_individual	-0.034	-0.109	0.041	0.376	-0.803	0.736
M_offspring	0.061	-0.063	0.185	0.333	-0.715	0.837
B_individual-B_offspring	0.098	-0.198	0.395	0.515	-0.723	0.920
B_individual-F_individual	0.067	-0.028	0.162	0.168	-0.705	0.839
B_individual-F_offspring	0.057	-0.067	0.182	0.365	-0.719	0.833
B_individual-M_individual	0.045	-0.053	0.142	0.370	-0.728	0.817
B_individual-M_offspring	0.140	0.002	0.277	0.047	-0.639	0.918
B_offspring-F_individual	-0.031	-0.329	0.266	0.836	-0.853	0.790
B_offspring-F_offspring	-0.041	-0.349	0.268	0.795	-0.867	0.785
B_offspring-M_individual	-0.054	-0.352	0.244	0.723	-0.876	0.768
B_offspring-M_offspring	0.041	-0.273	0.355	0.797	-0.787	0.869
F_individual-F_offspring	-0.010	-0.112	0.093	0.855	-0.782	0.763
F_individual-M_individual	-0.022	-0.087	0.043	0.499	-0.791	0.746
F_individual-M_offspring	0.073	-0.046	0.192	0.232	-0.703	0.848
F_offspring-M_individual	-0.013	-0.118	0.092	0.810	-0.786	0.760
F_offspring-M_offspring	0.082	-0.051	0.215	0.227	-0.695	0.860
M_individual-M_offspring	0.095	-0.025	0.215	0.121	-0.680	0.870

Sex x Fitness (higher level)

dat$sex_fitness_higher_level <- as.factor(paste(dat$sex, dat$fitness_higher_level ,sep = "_"))

mod_sex_fitness <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_fitness_higher_level - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex_fitness)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-276.9777   553.9554   575.9554   625.5846   576.3548   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1094  0.3308    681     no         effectID 
sigma^2.2  0.0121  0.1102    203     no          paperID 
sigma^2.3  0.0287  0.1693    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 673) = 387707548.1980, p-val < .0001

Test of Moderators (coefficients 1:8):
F(df1 = 8, df2 = 673) = 0.7301, p-val = 0.6649

Model Results:

                                            estimate      se     tval   df 
sex_fitness_higher_levelB_reproduction       -0.0228  0.0717  -0.3185  673 
sex_fitness_higher_levelB_survival           -0.0948  0.0514  -1.8463  673 
sex_fitness_higher_levelF_indirect fitness   -0.2305  0.2539  -0.9078  673 
sex_fitness_higher_levelF_reproduction       -0.0027  0.0388  -0.0696  673 
sex_fitness_higher_levelF_survival           -0.0238  0.0436  -0.5454  673 
sex_fitness_higher_levelM_indirect fitness   -0.2293  0.2070  -1.1076  673 
sex_fitness_higher_levelM_reproduction       -0.0273  0.0427  -0.6406  673 
sex_fitness_higher_levelM_survival            0.0002  0.0450   0.0051  673 
                                              pval    ci.lb   ci.ub    
sex_fitness_higher_levelB_reproduction      0.7502  -0.1636  0.1179    
sex_fitness_higher_levelB_survival          0.0653  -0.1957  0.0060  . 
sex_fitness_higher_levelF_indirect fitness  0.3643  -0.7289  0.2680    
sex_fitness_higher_levelF_reproduction      0.9445  -0.0788  0.0734    
sex_fitness_higher_levelF_survival          0.5856  -0.1093  0.0618    
sex_fitness_higher_levelM_indirect fitness  0.2684  -0.6358  0.1772    
sex_fitness_higher_levelM_reproduction      0.5220  -0.1111  0.0564    
sex_fitness_higher_levelM_survival          0.9960  -0.0882  0.0886    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_sex_fitness)

   R2_marginal R2_conditional 
   0.007397922    0.276989436 

orchard_plot(mod_sex_fitness, mod = "sex_fitness_higher_level", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S24. Add text.

tab_sex_fitness <- all_models(mod_sex_fitness,  mod = "sex_fitness_higher_level", type = "normal")

# saving tab_sex_fitness as RDS
saveRDS(tab_sex_fitness, here("Rdata", "tab_sex_fitness.RDS"))

Table S18. Add text.

tab_sex_fitness <-readRDS(here("Rdata", "tab_sex_fitness.RDS"))

tab_sex_fitness

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
B_reproduction	-0.023	-0.164	0.118	0.749	-0.800	0.754
B_survival	-0.095	-0.196	0.006	0.066	-0.865	0.676
F_indirect fitness	-0.231	-0.732	0.270	0.365	-1.144	0.682
F_reproduction	-0.004	-0.081	0.073	0.922	-0.771	0.764
F_survival	-0.026	-0.112	0.061	0.558	-0.794	0.743
M_indirect fitness	-0.230	-0.638	0.178	0.269	-1.096	0.636
M_reproduction	-0.028	-0.112	0.056	0.513	-0.796	0.740
M_survival	-0.001	-0.089	0.088	0.988	-0.769	0.768
B_reproduction-B_survival	-0.072	-0.231	0.088	0.378	-0.852	0.708
B_reproduction-F_indirect fitness	-0.208	-0.726	0.310	0.431	-1.131	0.715
B_reproduction-F_reproduction	0.019	-0.132	0.170	0.803	-0.759	0.798
B_reproduction-F_survival	-0.003	-0.159	0.153	0.972	-0.782	0.777
B_reproduction-M_indirect fitness	-0.207	-0.636	0.223	0.345	-1.083	0.669
B_reproduction-M_reproduction	-0.005	-0.159	0.149	0.949	-0.784	0.774
B_reproduction-M_survival	0.022	-0.135	0.180	0.781	-0.757	0.802
B_survival-F_indirect fitness	-0.136	-0.644	0.371	0.598	-1.053	0.781
B_survival-F_reproduction	0.091	-0.017	0.199	0.098	-0.680	0.862
B_survival-F_survival	0.069	-0.045	0.183	0.233	-0.703	0.841
B_survival-M_indirect fitness	-0.135	-0.552	0.282	0.525	-1.005	0.735
B_survival-M_reproduction	0.067	-0.046	0.180	0.247	-0.705	0.839
B_survival-M_survival	0.094	-0.023	0.211	0.114	-0.678	0.867
F_indirect fitness-F_reproduction	0.227	-0.275	0.729	0.375	-0.687	1.141
F_indirect fitness-F_survival	0.205	-0.297	0.708	0.423	-0.709	1.119
F_indirect fitness-M_indirect fitness	0.001	-0.602	0.604	0.997	-0.972	0.974
F_indirect fitness-M_reproduction	0.203	-0.300	0.706	0.428	-0.711	1.118
F_indirect fitness-M_survival	0.230	-0.272	0.733	0.369	-0.684	1.145
F_reproduction-F_survival	-0.022	-0.104	0.060	0.600	-0.790	0.746
F_reproduction-M_indirect fitness	-0.226	-0.636	0.184	0.279	-1.093	0.641
F_reproduction-M_reproduction	-0.024	-0.102	0.054	0.542	-0.792	0.743
F_reproduction-M_survival	0.003	-0.085	0.091	0.944	-0.766	0.772
F_survival-M_indirect fitness	-0.204	-0.615	0.207	0.330	-1.071	0.663
F_survival-M_reproduction	-0.002	-0.093	0.089	0.962	-0.771	0.767
F_survival-M_survival	0.025	-0.064	0.114	0.579	-0.744	0.794
M_indirect fitness-M_reproduction	0.202	-0.208	0.612	0.334	-0.665	1.069
M_indirect fitness-M_survival	0.229	-0.182	0.640	0.274	-0.638	1.096
M_reproduction-M_survival	0.027	-0.067	0.122	0.570	-0.742	0.797

Sex x Dispersal type

dat$sex_dispersal_type <- as.factor(paste(dat$sex, dat$dispersal_type ,sep = "_"))

mod_sex_dispersal_type <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_dispersal_type - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex_dispersal_type)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-277.1244   554.2489   578.2489   632.3720   578.7223   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1096  0.3311    681     no         effectID 
sigma^2.2  0.0157  0.1252    203     no          paperID 
sigma^2.3  0.0247  0.1570    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 672) = 387832987.1187, p-val < .0001

Test of Moderators (coefficients 1:9):
F(df1 = 9, df2 = 672) = 0.4659, p-val = 0.8977

Model Results:

                              estimate      se     tval   df    pval    ci.lb 
sex_dispersal_typeB_B          -0.0962  0.1144  -0.8409  672  0.4007  -0.3209 
sex_dispersal_typeB_breeding   -0.0040  0.0994  -0.0400  672  0.9681  -0.1991 
sex_dispersal_typeB_natal      -0.0834  0.0509  -1.6382  672  0.1018  -0.1835 
sex_dispersal_typeF_B          -0.0810  0.1115  -0.7259  672  0.4681  -0.2999 
sex_dispersal_typeF_breeding   -0.0032  0.0537  -0.0602  672  0.9520  -0.1086 
sex_dispersal_typeF_natal      -0.0061  0.0392  -0.1549  672  0.8770  -0.0830 
sex_dispersal_typeM_B          -0.1084  0.1228  -0.8826  672  0.3778  -0.3496 
sex_dispersal_typeM_breeding   -0.0215  0.0585  -0.3665  672  0.7141  -0.1364 
sex_dispersal_typeM_natal      -0.0030  0.0413  -0.0733  672  0.9416  -0.0842 
                               ci.ub    
sex_dispersal_typeB_B         0.1285    
sex_dispersal_typeB_breeding  0.1912    
sex_dispersal_typeB_natal     0.0166    
sex_dispersal_typeF_B         0.1380    
sex_dispersal_typeF_breeding  0.1022    
sex_dispersal_typeF_natal     0.0709    
sex_dispersal_typeM_B         0.1328    
sex_dispersal_typeM_breeding  0.0935    
sex_dispersal_typeM_natal     0.0781    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_sex_dispersal_type)

   R2_marginal R2_conditional 
   0.006633523    0.273905656 

orchard_plot(mod_sex_dispersal_type, mod = "sex_dispersal_type", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, trunk.size = 0.5, angle = 45)

Figure S25. Add text.

tab_sex_dispersal_type <- all_models(mod_sex_dispersal_type,  mod = "sex_dispersal_type", type = "normal")

# saving tab_sex_dispersal_type as RDS
saveRDS(tab_sex_dispersal_type, here("Rdata", "tab_sex_dispersal_type.RDS"))

Table S19. Add text.

tab_sex_dispersal_type <-readRDS(here("Rdata", "tab_sex_dispersal_type.RDS"))

tab_sex_dispersal_type

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
B_B	-0.096	-0.321	0.129	0.403	-0.892	0.699
B_breeding	-0.004	-0.200	0.192	0.967	-0.792	0.784
B_natal	-0.083	-0.184	0.017	0.103	-0.853	0.686
F_B	-0.081	-0.301	0.138	0.468	-0.875	0.713
F_breeding	-0.004	-0.110	0.102	0.940	-0.774	0.766
F_natal	-0.008	-0.086	0.070	0.840	-0.775	0.759
M_B	-0.108	-0.350	0.133	0.379	-0.909	0.692
M_breeding	-0.022	-0.137	0.094	0.712	-0.793	0.750
M_natal	-0.004	-0.086	0.077	0.922	-0.771	0.763
B_B-B_breeding	0.092	-0.203	0.387	0.540	-0.726	0.910
B_B-B_natal	0.013	-0.223	0.249	0.916	-0.786	0.811
B_B-F_B	0.015	-0.288	0.317	0.923	-0.806	0.836
B_B-F_breeding	0.092	-0.150	0.334	0.456	-0.708	0.892
B_B-F_natal	0.088	-0.142	0.318	0.453	-0.709	0.885
B_B-M_B	-0.012	-0.330	0.305	0.939	-0.839	0.814
B_B-M_breeding	0.074	-0.172	0.321	0.554	-0.727	0.876
B_B-M_natal	0.092	-0.139	0.323	0.435	-0.705	0.889
B_breeding-B_natal	-0.079	-0.292	0.133	0.464	-0.871	0.713
B_breeding-F_B	-0.077	-0.368	0.214	0.603	-0.893	0.739
B_breeding-F_breeding	0.000	-0.214	0.214	1.000	-0.792	0.792
B_breeding-F_natal	-0.004	-0.208	0.200	0.970	-0.794	0.786
B_breeding-M_B	-0.104	-0.413	0.204	0.506	-0.927	0.718
B_breeding-M_breeding	-0.018	-0.238	0.202	0.875	-0.812	0.776
B_breeding-M_natal	0.000	-0.206	0.206	1.000	-0.790	0.790
B_natal-F_B	0.002	-0.232	0.237	0.985	-0.796	0.800
B_natal-F_breeding	0.079	-0.051	0.209	0.232	-0.695	0.853
B_natal-F_natal	0.075	-0.032	0.183	0.169	-0.695	0.846
B_natal-M_B	-0.025	-0.281	0.231	0.848	-0.830	0.780
B_natal-M_breeding	0.062	-0.078	0.202	0.386	-0.714	0.837
B_natal-M_natal	0.079	-0.031	0.190	0.158	-0.692	0.850
F_B-F_breeding	0.077	-0.160	0.314	0.523	-0.722	0.876
F_B-F_natal	0.073	-0.153	0.300	0.526	-0.723	0.869
F_B-M_B	-0.027	-0.288	0.234	0.837	-0.834	0.779
F_B-M_breeding	0.059	-0.183	0.302	0.630	-0.741	0.860
F_B-M_natal	0.077	-0.150	0.304	0.506	-0.719	0.873
F_breeding-F_natal	-0.004	-0.113	0.105	0.944	-0.775	0.767
F_breeding-M_B	-0.104	-0.362	0.154	0.427	-0.910	0.701
F_breeding-M_breeding	-0.018	-0.130	0.095	0.758	-0.789	0.754
F_breeding-M_natal	0.000	-0.112	0.112	1.000	-0.771	0.771
F_natal-M_B	-0.100	-0.349	0.148	0.427	-0.903	0.702
F_natal-M_breeding	-0.014	-0.134	0.107	0.823	-0.786	0.759
F_natal-M_natal	0.004	-0.068	0.076	0.916	-0.762	0.770
M_B-M_breeding	0.087	-0.176	0.350	0.517	-0.720	0.894
M_B-M_natal	0.104	-0.145	0.354	0.412	-0.698	0.907
M_breeding-M_natal	0.018	-0.106	0.142	0.780	-0.755	0.791

5 Multi-moderator-regression model

Full model

# used "maximum likelihood" for model selection

#######################
# Multi-variable models
#######################

mod_full <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ 1 + sex +
                age_class_clean + whose_fitness + 
                fitness_higher_level +
                dispersal_type + dispersal_phase,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              method = "ML",
              test = "t",
              sparse = TRUE)
summary(mod_full)

Multivariate Meta-Analysis Model (k = 681; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-273.7918  9634.4147   575.5837   638.9136   576.2143   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1078  0.3283    681     no         effectID 
sigma^2.2  0.0076  0.0872    203     no          paperID 
sigma^2.3  0.0321  0.1793    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 670) = 382387194.5207, p-val < .0001

Test of Moderators (coefficients 2:11):
F(df1 = 10, df2 = 670) = 1.1948, p-val = 0.2909

Model Results:

                                  estimate      se     tval   df    pval 
intrcpt                            -0.2609  0.1981  -1.3168  670  0.1884 
sexF                                0.0460  0.0501   0.9191  670  0.3584 
sexM                                0.0427  0.0508   0.8405  670  0.4010 
age_class_cleanjuvenile             0.1029  0.0575   1.7899  670  0.0739 
age_class_cleanmix                 -0.0880  0.0680  -1.2941  670  0.1961 
whose_fitnessfocal                 -0.0516  0.0480  -1.0756  670  0.2825 
fitness_higher_levelreproduction    0.2118  0.1699   1.2467  670  0.2129 
fitness_higher_levelsurvival        0.1810  0.1725   1.0493  670  0.2944 
dispersal_typebreeding              0.0209  0.0844   0.2476  670  0.8045 
dispersal_typenatal                 0.0051  0.0771   0.0668  670  0.9468 
dispersal_phasesettlement           0.0379  0.0536   0.7066  670  0.4800 
                                    ci.lb   ci.ub    
intrcpt                           -0.6498  0.1281    
sexF                              -0.0523  0.1444    
sexM                              -0.0571  0.1425    
age_class_cleanjuvenile           -0.0100  0.2158  . 
age_class_cleanmix                -0.2216  0.0455    
whose_fitnessfocal                -0.1459  0.0426    
fitness_higher_levelreproduction  -0.1218  0.5454    
fitness_higher_levelsurvival      -0.1577  0.5197    
dispersal_typebreeding            -0.1448  0.1867    
dispersal_typenatal               -0.1463  0.1566    
dispersal_phasesettlement         -0.0673  0.1430    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

round(r2_ml(mod_full)*100, 2)

   R2_marginal R2_conditional 
          1.96          28.37 

Multi-model inference

# multi-model selection
candidates <- dredge(mod_full, trace = 2)

# saving candidates as RDS
saveRDS(candidates, here("Rdata", "candidates.RDS"))

candidates <-readRDS(here("Rdata", "candidates.RDS"))

candidates

Global model call: rma.mv(yi = yi, V = VCV, mods = ~1 + sex + age_class_clean + 
    whose_fitness + fitness_higher_level + dispersal_type + dispersal_phase, 
    random = list(~1 | effectID, ~1 | paperID, ~1 | species_cleaned), 
    data = dat, method = "ML", test = "t", sparse = TRUE)
---
Model selection table 
   (Int) age_cls_cln dsp_phs dsp_typ ftn_hgh_lvl sex whs_ftn df   logLik  AICc
2      +           +                                          6 -276.828 565.8
4      +           +       +                                  7 -275.933 566.0
34     +           +                                       +  7 -276.278 566.7
36     +           +       +                               +  8 -275.479 567.2
1      +                                                      4 -279.768 567.6
42     +           +                           +           +  9 -274.701 567.7
10     +           +                           +              8 -275.782 567.8
18     +           +                               +          8 -275.927 568.1
12     +           +       +                   +              9 -275.053 568.4
33     +                                                   +  5 -279.217 568.5
3      +                   +                                  5 -279.276 568.6
44     +           +       +                   +           + 10 -274.319 569.0
20     +           +       +                       +          9 -275.359 569.0
50     +           +                               +       +  9 -275.518 569.3
6      +           +               +                          8 -276.641 569.5
9      +                                       +              6 -278.813 569.8
35     +                   +                               +  6 -278.815 569.8
41     +                                       +           +  7 -277.874 569.9
8      +           +       +       +                          9 -275.830 569.9
17     +                                           +          6 -278.913 570.0
26     +           +                           +   +         10 -274.949 570.2
52     +           +       +                       +       + 10 -275.002 570.3
38     +           +               +                       +  9 -276.064 570.4
5      +                           +                          6 -279.306 570.7
58     +           +                           +   +       + 11 -274.187 570.8
40     +           +       +       +                       + 10 -275.349 571.0
11     +                   +                   +              7 -278.449 571.1
49     +                                           +       +  7 -278.515 571.2
28     +           +       +                   +   +         11 -274.474 571.3
46     +           +               +           +           + 11 -274.557 571.5
19     +                   +                       +          7 -278.735 571.6
14     +           +               +           +             10 -275.674 571.7
37     +                           +                       +  7 -278.774 571.7
43     +                   +                   +           +  8 -277.769 571.8
22     +           +               +               +         10 -275.770 571.9
7      +                   +       +                          7 -278.938 572.0
25     +                                       +   +          8 -278.038 572.3
60     +           +       +                   +   +       + 12 -273.929 572.3
16     +           +       +       +           +             11 -274.996 572.4
48     +           +       +       +           +           + 12 -274.220 572.9
24     +           +       +       +               +         11 -275.261 572.9
51     +                   +                       +       +  8 -278.367 572.9
54     +           +               +               +       + 11 -275.338 573.1
57     +                                       +   +       +  9 -277.439 573.1
13     +                           +           +              8 -278.487 573.2
39     +                   +       +                       +  8 -278.490 573.2
21     +                           +               +          8 -278.571 573.4
45     +                           +           +           +  9 -277.602 573.5
27     +                   +                   +   +          9 -277.898 574.1
30     +           +               +           +   +         12 -274.856 574.2
56     +           +       +       +               +       + 12 -274.881 574.2
53     +                           +               +       +  9 -278.186 574.6
15     +                   +       +           +              9 -278.197 574.7
62     +           +               +           +   +       + 13 -274.062 574.7
59     +                   +                   +   +       + 10 -277.406 575.1
23     +                   +       +               +          9 -278.444 575.2
32     +           +       +       +           +   +         13 -274.418 575.4
47     +                   +       +           +           + 10 -277.534 575.4
29     +                           +           +   +         10 -277.795 575.9
64     +           +       +       +           +   +       + 14 -273.837 576.3
55     +                   +       +               +       + 10 -278.088 576.5
61     +                           +           +   +       + 11 -277.225 576.8
31     +                   +       +           +   +         11 -277.687 577.8
63     +                   +       +           +   +       + 12 -277.207 578.9
   delta weight
2   0.00  0.122
4   0.25  0.107
34  0.94  0.076
36  1.39  0.061
1   1.81  0.049
42  1.89  0.047
10  2.00  0.045
18  2.29  0.039
12  2.59  0.033
33  2.74  0.031
3   2.86  0.029
44  3.19  0.025
20  3.21  0.024
50  3.52  0.021
6   3.72  0.019
9   3.97  0.017
35  3.97  0.017
41  4.13  0.015
8   4.15  0.015
17  4.17  0.015
26  4.45  0.013
52  4.55  0.012
38  4.62  0.012
5   4.96  0.010
58  4.99  0.010
40  5.25  0.009
11  5.28  0.009
49  5.42  0.008
28  5.56  0.008
46  5.73  0.007
19  5.86  0.007
14  5.90  0.006
37  5.93  0.006
43  5.97  0.006
22  6.09  0.006
7   6.26  0.005
25  6.51  0.005
60  6.55  0.005
16  6.61  0.004
48  7.13  0.003
24  7.14  0.003
51  7.17  0.003
54  7.29  0.003
57  7.37  0.003
13  7.41  0.003
39  7.42  0.003
21  7.58  0.003
45  7.69  0.003
27  8.28  0.002
30  8.40  0.002
56  8.45  0.002
53  8.86  0.001
15  8.88  0.001
62  8.89  0.001
59  9.36  0.001
23  9.38  0.001
32  9.60  0.001
47  9.62  0.001
29 10.14  0.001
64 10.53  0.001
55 10.72  0.001
61 11.07  0.000
31 11.99  0.000
63 13.10  0.000
Models ranked by AICc(x) 

# displays delta AICc <2
candidates_aic2 <- subset(candidates, delta < 2) 
# model averaging
mr_averaged_aic2 <- summary(model.avg(candidates, delta < 2)) 

mr_averaged_aic2

Call:
model.avg(object = candidates, subset = delta < 2)

Component model call: 
rma.mv(yi = yi, V = VCV, mods = ~<7 unique rhs>, random = list(~1 | 
     effectID, ~1 | paperID, ~1 | species_cleaned), data = dat, method = ML, 
     test = t, sparse = TRUE)

Component models: 
       df  logLik   AICc delta weight
1       6 -276.83 565.78  0.00   0.24
12      7 -275.93 566.03  0.25   0.21
14      7 -276.28 566.72  0.94   0.15
124     8 -275.48 567.17  1.39   0.12
(Null)  4 -279.77 567.59  1.81   0.10
134     9 -274.70 567.67  1.89   0.09
13      8 -275.78 567.78  2.00   0.09

Term codes: 
     age_class_clean      dispersal_phase fitness_higher_level 
                   1                    2                    3 
       whose_fitness 
                   4 

Model-averaged coefficients:  
(full average) 
                                 Estimate Std. Error z value Pr(>|z|)
intrcpt                          -0.08622    0.11029   0.782    0.434
age_class_cleanjuvenile           0.06241    0.05321   1.173    0.241
age_class_cleanmix               -0.10467    0.07067   1.481    0.139
dispersal_phasesettlement         0.02130    0.04137   0.515    0.607
whose_fitnessoffspring            0.01788    0.03542   0.505    0.614
fitness_higher_levelreproduction  0.03936    0.11059   0.356    0.722
fitness_higher_levelsurvival      0.03285    0.10117   0.325    0.745
 
(conditional average) 
                                 Estimate Std. Error z value Pr(>|z|)  
intrcpt                          -0.08622    0.11029   0.782   0.4344  
age_class_cleanjuvenile           0.06911    0.05169   1.337   0.1812  
age_class_cleanmix               -0.11591    0.06502   1.783   0.0746 .
dispersal_phasesettlement         0.06429    0.04902   1.312   0.1896  
whose_fitnessoffspring            0.04929    0.04369   1.128   0.2592  
fitness_higher_levelreproduction  0.21658    0.17002   1.274   0.2027  
fitness_higher_levelsurvival      0.18075    0.17198   1.051   0.2933  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# relative importance of each predictor for all the models
importance <- sw(candidates)

importance

                     age_class_clean dispersal_phase whose_fitness
Sum of weights:      0.74            0.40            0.40         
N containing models:   32              32              32         
                     fitness_higher_level sex  dispersal_type
Sum of weights:      0.28                 0.20 0.14          
N containing models:   32                   32   32          

6 Publication bias

Egger regression: uni-moderator

# if each group n is available - assume n/2
dat$effectN <- ifelse(is.na(dat$n_group_1), (dat$n/2)*2/dat$n,  
                      (dat$n_group_1 * dat$n_group_2) / (dat$n_group_1 + dat$n_group_2))
  
dat$sqeffectN <- sqrt(dat$effectN)

mod_egger <- rma.mv(yi = yi, 
                    V = VCV,
                mod = ~ sqeffectN,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                sparse = TRUE)
summary(mod_egger)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.7841   557.5682   567.5682   590.1714   567.6574   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1096  0.3311    681     no         effectID 
sigma^2.2  0.0117  0.1082    203     no          paperID 
sigma^2.3  0.0281  0.1677    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 679) = 390714012.2193, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 679) = 0.0041, p-val = 0.9487

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt     -0.0307  0.0403  -0.7624  679  0.4461  -0.1097  0.0484    
sqeffectN    0.0003  0.0045   0.0643  679  0.9487  -0.0085  0.0091    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

round(r2_ml(mod_egger)*100, 2)

   R2_marginal R2_conditional 
          0.00          26.65 

small <- bubble_plot(mod_egger,
                     mod = "sqeffectN",
                     group = "paperID",
                     xlab = "sqrt(Effective N)",
                     ylab = "Effect Size: Zr",
                     g = TRUE)

#| fig-cap: "**Figure S27.** Add text."

small

Decline effect: uni-moderator

# creating publication year from "reference"
dat$reference <- as.character(dat$reference)
dat$year <- with(dat, substr(reference, nchar(reference)-4, nchar(reference)))
dat$year <- as.integer(ifelse(dat$year == "2017a" | dat$year == "2017b", 2017, dat$year))
# decline effect
mod_dec <- rma.mv(yi = yi, V = vi,
                     mod = ~ year,
                     data = dat, 
                     random = list(
                       ~ 1 | effectID,
                       ~ 1 | paperID,
                       ~ 1 | species_cleaned),
                     test = "t",
                     sparse = TRUE)
summary(mod_dec)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-223.1394   446.2787   456.2787   478.8818   456.3679   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0630  0.2510    681     no         effectID 
sigma^2.2  0.0125  0.1116    203     no          paperID 
sigma^2.3  0.0268  0.1638    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 679) = 7844.2279, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 679) = 0.6237, p-val = 0.4299

Model Results:

         estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt   -2.4358  3.0575  -0.7967  679  0.4259  -8.4391  3.5675    
year       0.0012  0.0015   0.7898  679  0.4299  -0.0018  0.0042    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

round(r2_ml(mod_dec)*100, 2)

   R2_marginal R2_conditional 
          0.22          38.54 

decline <- bubble_plot(mod_dec,
                       mod = "year",
                       group = "paperID",
                       xlab = "Publication year",
                       ylab = "Effect Size: Zr",
                       g = TRUE)

decline 

Figure S28. Add text.

Combined analysis

# All together

mod_comb <- rma.mv(yi = yi, 
                   V = VCV,
                   mod = ~ year + sqeffectN,
                   data = dat, 
                   random = list(
                     ~ 1 | effectID,
                     ~ 1 | paperID,
                     ~ 1 | species_cleaned),
                   # ~ 1 | phylogeny),
                   # R= list(phylogeny = cor_tree),
                   test = "t",
                   sparse = TRUE)
                   
summary(mod_comb)

Multivariate Meta-Analysis Model (k = 681; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.3077   556.6154   568.6154   595.7303   568.7406   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1092  0.3304    681     no         effectID 
sigma^2.2  0.0110  0.1046    203     no          paperID 
sigma^2.3  0.0308  0.1754    147     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 678) = 390510267.1254, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 678) = 0.1988, p-val = 0.8198

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt     -2.1383  3.3535  -0.6376  678  0.5239  -8.7229  4.4462    
year         0.0011  0.0017   0.6282  678  0.5301  -0.0022  0.0043    
sqeffectN   -0.0002  0.0045  -0.0478  678  0.9619  -0.0091  0.0087    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

round(r2_ml(mod_comb)*100, 2)

   R2_marginal R2_conditional 
          0.11          27.72 

# small-study
bubble_plot(mod_egger,
            mod = "sqeffectN",
            group = "paperID",
            xlab = "sqrt(Effective N)",
            ylab = "Effect Size: Zr",
            g = TRUE)

Figure S29. Add text.

# decline
bubble_plot(mod_comb,
            mod = "year",
            group = "paperID",
            xlab = "Publication year",
            ylab = "Effect Size: Zr",
            g = TRUE)

Figure S30. Add text.

Leave-1study-out (sensitivity analysis)

dat <- dat %>%
  mutate(leave_out = reference)

dat$leave_out <- as.factor(dat$leave_out)


LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$leave_out))) {
  temp_dat <- dat %>%
    filter(leave_out != levels(dat$leave_out)[i])
  
  VCV_leaveout <- vcalc(vi = temp_dat$vi, cluster = temp_dat$shared_group, rho = 0.5)
  
  LeaveOneOut_effectsize[[i]] <-  rma.mv(yi = yi,
                                         V = VCV_leaveout, 
                                         random = list(
                                           ~ 1 | effectID,
                                           ~ 1 | paperID,
                                           ~ 1 | species_cleaned),
                                         test = "t",
                                         sparse = TRUE,
                                         data = temp_dat[temp_dat$leave_out != levels(temp_dat$leave_out)[i], ])
}

# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(model) {
  df <- data.frame(est = model$b, lower = model$ci.lb, upper = model$ci.ub)
  return(df)
}

# using dplyr to form data frame
MA_oneout <- lapply(LeaveOneOut_effectsize,function(x) est.func(x)) %>%
  bind_rows %>%
  mutate(left_out = levels(dat$leave_out))


# telling ggplot to stop reordering factors
MA_oneout$left_out <- as.factor(MA_oneout$left_out)
MA_oneout$left_out <- factor(MA_oneout$left_out, levels = MA_oneout$left_out)

# saving the runs
saveRDS(MA_oneout, here("Rdata", "MA_oneout.RDS"))

MA_oneout <- readRDS(here("Rdata", "MA_oneout.RDS"))

# plotting
leaveoneout <- ggplot(MA_oneout) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +

  geom_hline(yintercept = mod1$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
  geom_hline(yintercept = mod1$b, lty = 1, lwd = 0.75, colour = "black") + 
  geom_hline(yintercept = mod1$ci.ub,
             lty = 3, lwd = 0.75, colour = "black") + 
  geom_pointrange(aes(x = left_out, y = est,
                      ymin = lower, ymax = upper)) + 
  xlab("Study left out") + 
  ylab("Zr (effect size), 95% CI") +
  coord_flip() + 
  theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
  theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))

leaveoneout

Figure S31. Add text.

7 Figure gallery

Figure 1

reduced_intercept

Test caption

Figure —

8 R Session Information

# pander for making it look nicer
sessionInfo() %>% pander()

R version 4.2.0 (2022-04-22)

Platform: x86_64-apple-darwin17.0 (64-bit)

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: matrixcalc(v.1.0-6), orchaRd(v.2.0), here(v.1.0.1), patchwork(v.1.1.3), ape(v.5.7-1), metafor(v.4.4-0), numDeriv(v.2016.8-1.1), metadat(v.1.2-0), Matrix(v.1.5-4.1), pander(v.0.6.5), gridExtra(v.2.3), emmeans(v.1.8.9), MuMIn(v.1.47.5), kableExtra(v.1.3.4), lubridate(v.1.9.3), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.3), purrr(v.1.0.2), readr(v.2.1.4), tidyr(v.1.3.0), tibble(v.3.2.1), ggplot2(v.3.4.4) and tidyverse(v.2.0.0)

loaded via a namespace (and not attached): httr(v.1.4.7), jsonlite(v.1.8.7), viridisLite(v.0.4.2), splines(v.4.2.0), highr(v.0.10), stats4(v.4.2.0), vipor(v.0.4.7), yaml(v.2.3.7), pillar(v.1.9.0), lattice(v.0.22-5), glue(v.1.6.2), digest(v.0.6.33), rvest(v.1.0.3), colorspace(v.2.1-0), sandwich(v.3.1-0), htmltools(v.0.5.6.1), pkgconfig(v.2.0.3), xtable(v.1.8-4), mvtnorm(v.1.2-3), scales(v.1.3.0), webshot(v.0.5.5), svglite(v.2.1.2), tzdb(v.0.4.0), timechange(v.0.2.0), mgcv(v.1.9-0), farver(v.2.1.1), generics(v.0.1.3), TH.data(v.1.1-2), pacman(v.0.5.1), withr(v.2.5.2), cli(v.3.6.1), survival(v.3.5-7), magrittr(v.2.0.3), estimability(v.1.4.1), evaluate(v.0.23), fansi(v.1.0.5), nlme(v.3.1-163), MASS(v.7.3-60), xml2(v.1.3.5), beeswarm(v.0.4.0), tools(v.4.2.0), hms(v.1.1.3), lifecycle(v.1.0.4), multcomp(v.1.4-25), munsell(v.0.5.0), compiler(v.4.2.0), systemfonts(v.1.0.5), rlang(v.1.1.1), grid(v.4.2.0), rstudioapi(v.0.15.0), htmlwidgets(v.1.6.1), labeling(v.0.4.3), rmarkdown(v.2.25), gtable(v.0.3.4), codetools(v.0.2-19), R6(v.2.5.1), zoo(v.1.8-12), knitr(v.1.45), fastmap(v.1.1.1), utf8(v.1.2.3), rprojroot(v.2.0.4), mathjaxr(v.1.6-0), latex2exp(v.0.9.6), ggbeeswarm(v.0.7.2), stringi(v.1.7.12), parallel(v.4.2.0), Rcpp(v.1.0.11), vctrs(v.0.6.4), tidyselect(v.1.2.0), xfun(v.0.40) and coda(v.0.19-4)